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MILITARY  HYDROLOGY 


BREACH  EROSION  OF  EARTH-FILL  DAMS  AND 
FLOOD  ROUTING  (BEED)  MODEL 


PART  I :  INTRODUCTION 


Background 

1.  Under  the  Meteorological/Environmental  Plan  for  Action,  Phase  II, 
approved  for  implementation  on  26  January  1983,  the  US  Army  Corps  of  Engi¬ 
neers  (USACE)  has  been  tasked  to  implement  a  Research,  Development,  Testing, 
and  Evaluation  program  that  will  (a)  provide  the  Army  with  environmental 
effects  information  needed  to  operate  in  a  realistic  battlefield  environment 
and  (b)  provide  the  Army  with  the  capability  for  near-real  time  environmental 
effects  assessment  on  military  materiel  and  operations  in  combat.  In  response 
to  this  tasking,  the  Directorate  for  Research  and  Development,  USACE,  initi¬ 
ated  the  AirLand  Battlefield  Environment  (ALBE)  Thrust  program.  This  new  ini¬ 
tiative  will  develop  the  technologies  to  provide  the  field  Army  with  the 
operational  capability  to  perform  and  exploit  battlefield  effects  assessments 
for  tactical  advantage. 

2.  Military  hydrology,  one  facet  of  the  ALBE  Thrust,  is  a  specialized 
field  ot  study  tnat  deals  with  the  efleccs  o I  suriace  and  subsurface  water  on 
planning  and  conducting  military  operations.  In  1977,  the  Headquarters, 

USACE,  approved  a  military  hydrology  research  program;  management  responsi¬ 
bility  was  subsequently  assigned  to  the  Environmental  Laboratory,  US  Army 
Engineer  Waterways  Experiment  Station,  Vicksburg,  MS. 

3.  The  objective  of  military  hydrology  research  is  to  develop  an 
improved  hydrologic  capability  for  the  Armed  Forces  with  emphasis  on  applica¬ 
tions  in  the  tactical  environment.  To  meet  this  overall  objective,  research 
is  being  conducted  in  four  areas:  (a)  weather-hydrology  interactions, 

(b)  state  of  the  ground,  (c)  streamflow,  and  (d)  water  supply. 

4.  Previously  published  Military  Hydrology  reports  are  lifted  inside  the 
back  cover.  This  report  is  the  fifth  that  contributes  to  the  streamflow 
modeling  area.  Streamflow  modeling  is  oriented  toward  the  development  of  pro¬ 
cedures  for  rapidly  forecasting  streamflow  parameters,  including  discharge, 


velocity,  depth,  width,  and  flooded  area  from  natural  and  man-induced 
hydrologic  events.  Specific  work  efforts  include  (a)  the  development  of 
simple  and  objective  streamflow  forecasting  procedures  suitable  for  Armv  Ter¬ 
rain  Team  use,  (b)  the  adaptation  of  procedures  to  automatic  data  processing 
equipment  available  to  Terrain  Teams,  (c)  the  development  of  procedures  for 
accessing  and  processing  information  included  in  digital  terrain  data  bases, 
and  (d)  the  development  of  streamflow  analysis  and  display  concepts. 

Purpose  and  Scope 

5.  The  worn  reported  herein  is  an  effort  in  the  "Induced  Floods  as 
Linear/Area  Obstacles"  work  unit  of  Department  of  the  Armv  Project 

No.  4A7627 1 9AT40 .  The  objective  of  the  work  unit  is  to  provide  the  Armed 
Forces  improved  capabilities  for  forecasting  the  downstream  flood  flow  impacts 
resulting  from  controlled  or  uncontrolled  (dam  breach)  releases  for  single  or 
multiple  dams. 

6.  The  purpose  of  this  study  was  threefold: 

a.  To  develop  a  mathematical  model  for  the  simulation  of  gradual 
erosion  processes  of  an  earth  dam  so  that  the  flash-flood  hydro¬ 
graph  can  he  predicted. 

b.  To  route  the  released  water  mass  through  a  certain  distance  down¬ 
stream  by  means  of  an  existing  numerical  technique. 

c.  To  conduct  a  sensitivity  analysis  for  the  various  parameters 
involved . 

7.  The  first  phase  included  development  of  a  numerical  model,  both  for 
mainframe  and  microcomputer  facilities,  as  well  as  analytical  solutions  for 
simplified  versions  thereof  for  prediction  of  the  flash-flood  hvdrograph.  In 
the  second  phase,  the  solutio--  provided  in  the  first  phase  were  used  as 
upstream  boundary  conditions  for  the  Musk ingum-Cunge  method  with  variable 
parameters  that  will  route  the  flood  wave  through  the  receiving  downstream 
channel.  In  the  third  phase,  the  combined  model  was  applied  under  various 
conditions,  and  the  results  were  compared  and  analyzed. 

Dam  Failure 

8.  Devastating  flash  floods  resulting  t rom  sudden  dam  failure  involve 
potential  hazard  to  both  human  life  and  property.  Jansen  (1980)  states  that 
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there  havi-  jeen  approximately  2,000  dam  failures  around  the  world  since  the 
12th  .tury.  About  10  percent  of  those  failures  occurred  during  the  20th 
century,  causing  loss  of  more  than  8  000  lives  and  damage  costs  of  millions  of 
dollars.  A  recent  example  is  the  failure  of  Stava  Valley  Dam  in  Italy  on 
10  duly  1985,  which  resulted  in  200  fatalities  and  the  destruction  of  20 
houses  and  3  hotels. 

9.  The  International  Commission  on  Large  Dams  census  of  1962  registered 

9,315  dams  with  heights  greater  than  15  m  or  between  10  and  15  m  if  water 
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storage  exceeds  1,000,000  m  .  However,  according  to  Gruner  (1967),  the  total 
number  of  dams  that  impose  risk  of  serious  damage  in  case  of  failure  may  well 
exceed  150,000.  The  I'S  Army  Corps  of  Engineers  (1975)  classified  about  20,000 
dams  in  the  I'nited  States  as  potentially  dangerous  in  the  event  of  a  failure. 
In  spite  of  these  impressive  statistics,  little  is  known  about  the  triggering 
and  controlling  mechanisms  of  dam  failure. 

10.  The  majority  of  dams  are  man-made  earth-filled  dams.  Their  failure 
can  be  attributed  to  a  single  factor  or  to  a  combination  of  various  factors 
such  as  unexpectedly  large  inflows,  inadequate  foundation,  differential 
settlement,  landslides,  earthquakes,  poor  design  or  r onst rue t ion  ,  deficient 
materials,  improper  management,  or  acts  of  war.  The  mode  of  failure  depends 
both  on  the  cause  and  rhe  characteristics  of  the  individual  dam.  Historical 
earth-dam  failure  data  indicate  that  the  time  taken  for  the  reservoir  to  empty 
after  the  dam  was  breached  has  varied  from  15  min  to  more  than  5  hr  (Singh  and 
Snorrason  19821.  This  is  an  indication  that  dam  failure  is  a  time-dependent 
and  not  an  instantaneous  process. 

11.  The  failure  processes  on  an  earth  dam  are  generally  classified  in  one 

of  the  following  categories:  internal  erosion  due  to  piping,  progressive  ero¬ 

sion  of  the  downstream  face  due  to  seepage,  or  overtopping  of  the  crest  and 
subsequent  enlargement  from  erosion  of  an  initially  developed  breach.  Statis¬ 
tics  based  on  information  from  several  sources  (l.ou  1981)  show  that  about 

40  percent  of  failures  are  caused  bv  piping  or  seepage,  30  percent  by  overtop¬ 
ping  due  to  inadequate  spillway  capacity,  10  nercent  bv  landslides,  and  20 
percent  bv  other  causes  classified  as  miscellaneous.  The  ability  to  predict 
dam  bleaching  is  essential  for  a  reliable  estimation  of  the  released  water 
hvdrograph.  The  shape,  duration,  and  magnitude  of  the  d^m- breach  flood  hydro¬ 
graph  affect  the  results  of  flood  routing  on  its  downstream  course.  Accuracy 


of  these  results  is  very  important  for  flood  forecasting,  contingency  evacua¬ 
tion  planning,  and  management  decisions  for  dam  safety. 

12.  The  dam-hreak  problem  can  be  divided  into  two  parts:  dam  failure 
processes  and  routing  of  the  released  mass  of  water  downstream.  l'he  two  pacts 
can  be  solved  separately.  Of  course,  rhe  sequence  of  the  solution  cannot  be 
changed,  since  the  results  of  the  first  part  must  be  used  as  upstream  boundary 
conditions  for  the  study  of  the  second  part. 

13.  Failure  of  an  earth  dam  is  a  very  complicated,  unsteady,  nonhomoge- 
neous,  three-dimensional  phenomenon  that  is  still  not  fully  understood.  Ihe 
size,  shape,  and  location  of  the  initial  breach  is  usually  unknown.  The  ero¬ 
sion  processes  during  breach  enlargement  involving  suspended  sediment  trans¬ 
port,  layer  by  layer  and/or  mass  erosion,  and  sloughing  of  the  slopes  are  verv 
dynamic  processes  that  have  not  been  defined  theoretically  as  yet.  On  the 
other  hand,  routing  of  the  flood  wave  downstream  becomes  complicated  bv  rapid 
changes  of  the  morphology  of  the  receiving  channel  or  ’  asin  due  to  scour' ng  or 
shoaling,  inadequate  information  regarding  friction  factors,  and  water  mass 
losses  due  to  infiltration  or  local  storage.  Another  process  contributing  to 
the  complexity  of  the  problem  is  the  presence  of  solid  materials,  in  the  term 
of  mud  and  debris,  which  are  carried  downstream  by  the  flowing  water. 

14.  In  spite  of  these  difficulties,  it  is  possible  to  idealize  the  system 
and  to  develop  a  mathematical  model  for  dam  break/flood  routing  simulation  by 
making  proper  assumptions  and  simplifications.  The  accuracy  of  this  model 
will  be  compatible  with  the  validity  of  its  approximations.  Due  to  the  large 
number  of  controlling  physical  parameters,  the  uncertainty  of  the  governing 
processes,  and  the  idealization  of  the  physical  system,  it  is  essential  for 
the  sake  of  safety  to  predict  the  most  severe  conditions  to  be  expected  by 
conducting  a  sensitivity  analysis.  This  will  also  provide  information  about 
the  importance  of  each  individual  quantity  or  process  within  the  entire  dam 
break/flood  routing  simulation  model. 
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PART  II:  LITERATURE  REVIEW 


15.  Although  dam  break  and  flood  routing  might  be  coupled  processes,  they 
will  be  treated  separately  for  simplicity  throughout  this  study. 

Dam-Break  Mathematical  Modeling 


16.  In  spite  of  the  importance  of  the  subject,  very  few  attempts  have 
been  made  to  mathematically  model  the  gradual  failure  of  an  earth  dam.  All  of 
the  existing  models  are  based  on  the  princit>les  of  hvdraulics,  hydrodynamics, 
and  sediment  transport,  but  each  model  has  i:s  own  characteristic  features.  A 
general  discussion  of  these  models  is  giver  in  the  following  paragraphs. 
Cristof ano 

17.  The  first  attempt  to  simulate  the  mechanics  of  gradual  dam  breach 
erosion  was  perhaps  done  by  Cristofano  (1965) ,  who  equated  the  force  of  water 
flowing  over  the  breach  to  the  friction  resistance  force  acting  on  the  wetted 
perimeter  of  the  breach.  Alter  some  manipulation,  he  derived  a  differential 
equation  relating  the  rate  of  change  of  water  discharge  to  the  rate  of  change 
of  the  vertical  and  lateral  erosion  within  the  notch.  However,  the  applica¬ 
tion  of  this  equation  was  cumbersome  for  mrnual  computation  and  was  also 
discontinuous  in  certain  cases.  Cristofano  simplified  his  approach,  and  the 
following  analytical  expression  was  obtained: 
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water  discharge  through  the  breach 
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length  of  the  breach  in  the  flow  direction 
developed  angle  of  repose  of  the  soil 
hydraulic  head  at  any  given  time 
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18.  The  applicability  of  Cristofano’s  model  is  limited  by  the  assumption 

of  a  trapezoidal  breach  of  constant  width  where  the  side  slopes  and  the  longi¬ 
tudinal  slope  of  the  breach  bottom  are  equal  to  the  developed  angle  of 
internal  friction.  There  also  is  an  uncertainty  in  the  estimation  of  the  pro¬ 
portionality  constant  .  The  solution  requires  a  trial-and- error 

procedure.  The  model  was  applied  by  the  Bureau  of  Reclamation  to  Hyrum  Dam, 
Utah,  and  by  the  Tennessee  Valley  Authority  to  Brown’s  Ferry  Nuclear  Power 
Plant . 

Harris  and  Wagner 

19.  Harris  and  Wagner  (1967)  treated  the  dam  failure  problem  as  a  para¬ 
bolic  breach  subjected  to  erosion.  The  sediment  transport  was  estimated  by 
the  Schoklitsch  bed-load  formula.  The  flow  through  the  breach  was  assumed  as 
spillway  overflow,  while  tailwater  effects  were  neglected.  The  model  requires 
specification  of  breach  dimensions  and  slope,  in  addition  to  sediment  grain 
size  and  critical  value  of  discharge  for  Initiation  of  sediment  motion.  The 
applicability  of  the  model  is  limited  by  the  uncertainty  of  the  values  of  var¬ 
ious  parameters  involved  and  by  neglecting  tailwater  effects  and  sloughing. 
Brown  and  Rogers 

20.  Brown  and  Rogers  (1977,  1981)  reported  on  the  Bureau  of  Reclamation 
computer  program  BRDAM,  which  was  based  on  the  work  of  Harris  and  Wagner 
(1967).  The  model,  which  is  capable  of  simulating  erosion  from  either 
overtopping  or  piping,  was  applied  to  the  failure  of  Teton  Dam,  Idaho.  Its 
limitations  are  similar  to  those  of  the  original  model  of  Harris  and  Wagner. 
Fread 

21.  Extensive  research  on  dam-breach  flash  flooding  was  accomplished  by 
Fread  (1977,  1978,  1980,  1981)  using  the  National  Weather  Service  computer 
program  DAMBRK,  which  can  handle  rectangular,  triangular,  or  trapezoidal 
breach  shapes.  The  breaching  of  the  dam  commences  after  the  water  elevation 
within  the  reservoir  exceeds  a  specific  value,  and  the  breach  bottom  enlarges 
at  a  predetermined  linear  rate.  In  the  total  outflow  discharge,  both  broad- 
crested  weir  flow  over  the  breach  and  flow  through  spillway  outlets  are  incor¬ 
porated.  The  DAMBRK  program  was  applied  to  five  historic  dam-break  flood 
cases.  Although  the  results  after  calibration  were  satisfactory,  the  model 
cannot  be  applied  for  predictive  purposes  due  to  the  requirement  of  a  priori 
definitions  of  failure  time  duration  and  terminal  shape  and  size  of  the 
breach.  Therefore,  this  model  is  useful  only  for  the  estimation  of  a  spectrum 
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of  possible  flooding  events,  not  for  prediction  of  the  one  most  likely  to 


occur . 

Lou 

22.  Lou  (1981)  presented  a  model  for  estimation  of  the  outflow  hydrograph 
generated  by  a  gradual  earth-dam  rupture.  His  model  was  based  on  the  conti¬ 
nuity  and  momentum  equations  of  unsteady  flow  solved  by  Priessmann's  four- 
point  finite-difference  scheme.  The  inertia  terms  of  the  momentum  equation 
were  neglected.  For  the  sedimentation  processes,  Lou  initially  used  DuBoy's 
bed-load  equation  along  with  Einstein's  theory  for  suspended  sediment  trans¬ 
port.  However,  this  approach,  when  applied  to  dam-erosion  cases,  experienced 
instability  problems.  Thus,  he  proceeded  with  a  simplified  sediment  transport 
expression  that  he  called  a  transport  function.  It  was  derived  from  the 
assumption  that  embankment  erosion  was  proportional  to  the  kinetic  energy  of 
the  flowing  water  and  was  expressed  by  the  following  equation: 

M  =  eitdu4  (2) 


where 

M  =  mass  of  eroded  soil 

e.  =  erodibility  index 
x  J 

t,  =  failure  duration  time 
d 

u  =  water  velocity  through  the  breach 

23.  Applicability  of  the  model  as  a  predicter  is  very  limited  since  the 
duration  of  failure  time  and  the  erodibility  index  are  almost  impossible  to 
predetermine.  The  model  was  calibrated  and  tested  using  the  transport  func¬ 
tion  approach  for  the  failure  cases  of  Teton  Dam,  Idaho,  and  Mantaro  Dam, 

Peru.  The  results  were  satisfactory. 

Ponce  and  Tsivoglou 

24.  Ponce  and  Tsivoglou  (1981)  developed  a  gradual  dam-breach  model  using 
the  St.  Venant  system  of  equations,  which  they  solved  numerically  by  the 
Priessmann's  finite-difference  scheme.  The  sediment  routing  was  done  by  an 
Exner-type  equation  where  the  bed-load  function  was  that  of  Mever-Peter  and 
Mueller  (Simons  and  Senturk  1976).  Regarding  the  breach  morphology,  they 
introduced  a  regime-type  relation  between  top  width  of  the  breach  and  flow 
rate.  This  relation  was  applied  from  inception  to  peak  flow,  after  which  the 
breach  was  kept  constant.  The  weakness  of  this  model  is  the  determination  of 


the  rate  of  growth  of  the  breach  width  and  the  neglect  of  the  sloughing 
effects.  The  model  was  tested  on  actual  data  of  the  failure  of  the  natural 
embankment  that  formed  Mantaro  Dam  in  Peru. 


Fread 

25.  The  latest  development  on  breach  erosion  for  earth-fill  dams  is  the 
BREACH  model  presented  by  Fread  (1984).  This  is  an  iterative  numerical  model 
based  on  broad-crested  weir  flow  over  the  breach  and  quasi  steady-state  uni¬ 
form  flow  along  the  downstream  face  breach  channel.  In  development  of  the 
model,  tailwater  effects  were  included.  Sediment  transport  was  treated  by  the 
Meyer-Peter  and  Mueller  bed-load  formula.  The  innovative  aspect  of  the  model 
is  the  introduction  of  slope  stability,  although  the  theoretical  derivation  is 
for  dry  soil  conditions.  The  simulation  of  erosion  assumes  that  the  breach 
slope  is  parallel  to  the  downstream  face  slope  of  the  dam.  The  applicability 
of  the  model  for  predictive  purposes  is  restricted  by  the  uncertainty  of  the 
values  of  critical  shear  stress  for  initiation  of  erosion  and  terminal  breach 
width,  which  are  required  as  input  data  by  the  model.  The  model  was  applied 
to  the  failures  of  Teton  Dam,  Idaho,  and  Mantaro  Dam,  Peru. 

Classification  and  comparison  of  models 

26.  All  of  the  existing  models  have  some  advantages  and  disadvantages 
regarding  computational  efficiency  and  realistic  description  of  the  physical 
processes.  When  they  were  applied  to  historical  dam-failure  cases,  all  of  the 
models  showed  an  acceptable  degree  of  accuracy.  Of  course  this  is  due  par¬ 
tially  to  the  fact  that  a  number  of  parameters  can  be  calibrated  to  improve 
simulation  results.  The  basic  philosophy  for  mathematical  modeling  of  dam- 
break  problems  is  the  coupled  treatment  of  the  two  phases  involved,  i.e., 
reservoir  water  and  sediment  from  the  dam  body.  Governing  equations  and  the 
number  and  nature  of  physical  and  empirical  parameters  determine  the  suit¬ 
ability  of  the  model  for  prediction. 

27.  Similarities  and  differences  of  the  various  models  are  given  in 
Table  1.  This  illustrates  the  evolution  and  expansion  of  the  technology  of 
earth-dam  failure  simulation  during  the  last  20  years,  from  the  simple  con¬ 
ceptual  model  of  Cristofano  to  the  most  sophisticated  BREACH  model  of  Fread. 
Improvement  of  the  existing  state  of  the  art  can  be  achieved  by  reducing  the 
number  of  parameters  needing  calibration  and  by  introducing  more  realistic 
assumptions  for  both  water  discharge  and  sediment  transport  mechanisms. 


Mathemat ^ cal  Models  for  Dam-Breach  Erosion 
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Note.  Brae cet  indicates  that  models  are  included  in  the  same  category  due  to  their  simllariti 


Breach  Characteristics 


28.  One  of  the  weak  points  in  the  studies  of  earth-dam  failures  is  the 
breach  morphology.  Breach  shapes  and  dimensions  have  been  documented  in  many 
cases,  but  predictive  correlations  are  very  limited.  The  same  is  true 

of  the  failure  duration  time.  Information  on  pertinent  earth-dam  breach 
characteristics  for  52  cases  is  given  in  Table  2  (Ponce  1982;  Singh  and 
Snorrason  1982;  MacDonald  and  Langridge-Monopolis  1984). 

29.  Singh  and  Snorrason  (1982)  analyzed  the  historic  data  for  20  dams  and 
provided  information  on  the  three  breach  parameters:  width  of  breach,  initial 
hydraulic  head  for  failures  caused  by  overtopping,  and  failure  duration  time. 

30.  Ponce  (1982)  presented  a  preliminary  analysis  of  certain  parameters 
relevant  to  the  breach  morphology.  For  his  analysis  he  used  the  breach  Froude 
number  F  : 


F  = 


and  a  shape  factor  S^,  ,  defined  as 


S 


F 


Bd 

B  Z 
D  o 


(3) 


(4) 


where 

Q  =  peak  outflow  discharge 
B  =  top  width  of  the  breach 
g  =  acceleration  due  to  gravity 
d  =  depth  of  breach 

8^  =  top  width  of  dam 

Z  =  initial  height  of  dam 
o 

By  plotting  the  data  from  29  historical  cases  (Figure  1),  Ponce  derived  the 
relation 
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•  <re  presented  as  ap.s C  ream/downst  ream. 
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Figure  1.  Breach  Froude  number  versus  shape  factor 
(after  Ponce  1982) 

which  is  comparable  to  the  equation  reported  by  the  US  Army  Corps  of  Engineers 
(1961), 


F 


0.29S 


-0.28 

F 


(6) 


31.  Another  interesting  compilation  of  breach  characteristics  data  was 
presented  by  MacDonald  and  Langridge-Monopolis  (1984),  who  analyzed  42  cases 
and  suggested  an  empirical  methodology  for  predicting  the  shape,  size,  and 
failure  time  for  an  earth-fill  dam.  Their  methodology  is  based  on  Fig¬ 
ures  2-4.  In  Figures  2  and  3,  they  make  use  of  a  "breach  formation  factor," 
which  is  defined  as  the  product  of  the  discharged  volume  of  water  and  the 
difference  in  elevation  between  peak  reservoir  water  surface  and  breach  base. 
By  estimating  the  breach  formation  factor,  they  obtain  breach  volume  from 
Figure  2.  Having  the  volume  of  breach,  they  use  Figure  4  for  prediction  of 
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BREACH  FORMATION  FACTOR  IN  CU  METERS  x  METERS 


Figure  2.  Outflow  characteristics  versus  breach  size  (reference 
numbers  are  identified  in  Table  2)  (after  MacDonald  and  Langridge- 

Monopolis  1984) 
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BREACH  FORMATION  FACTOR  IN  CU  METERS  x  METERS 


VOLUME  OF  ERODED  DAM  MATERIALS  IN  CU  METERS 

Figure  4.  Breach  failure  time  versus  volume  (reference 
numbers  are  identified  in  Table  2)  (after  MacDonald  and 
Langr idge-Monopolis  1984) 

failure  time.  The  same  authors  suggested  a  triangular  breach  shape  with  2V:  1H 
side  slopes,  which  turns  into  a  trapezoidal  shape  after  the  breach  reaches  the 
base  of  the  dam.  Houston  (1984)  reanalyzed  the  previous  data,  proposing  a 
trapezoidal  breach  with  IV: 1H  side  slopes  and  base  width  equal  to  the  depth  of 
the  breach. 

32.  Further  analysis  of  breach  characteristics  is  given  in  Figures  5  and 
6  where  the  breach  top  (B) ,  bottom  (b) ,  and  average  widths  are  plotted  versus 
the  height  of  dam  and  the  depth  of  breach,  respectively.  Using  least  squares 
curve  fitting  approximation,  the  following  relations  were  obtained: 


B  =  4.43Z 

c 

b  =  2.C4Z 

c 

b  =  2 . 6Cd 
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BREACH  IN  METERS 


The  constant  coefficient  in  Equation  8  is  less  than  the  one  in  Equation  9 
because  breach  depth  d  is  sometimes  less  than  the  dam  height  ,  which 

means  partial  failure  occurred.  Indeed,  in  9  of  39  documented  dam  cases,  the 
failure  was  partial. 

33.  Based  m  the  data  of  Table  2,  the  probability  of  exceedance  of  dam 
failure  time  is  plotted  in  Figure  7.  In  the  same  figure,  the  probability  of 
exceedance  of  the  initial  hydraulic  head  for  an  overtopping  failure  event  is 
also  plotted  using  data  from  Singh  and  Snorrason  (1982).  Thus,  with  a 
50-peicent  probability,  failure  time  will  be  about  1.10  hr,  while  initial  head 
will  be  approxiately  0.4  m. 

34.  Although  these  results  provide  valuable  information  about  the  order 
of  magnitude  of  breach  characteristics  when  applied,  they  should  be  used  with 
caution  and  judgment.  The  scattering  of  the  data  points  and  lack  of 

INITIAL  OVERTOPPING  HEAD,  METERS 


Probability  of  exceedance  of  initial  overtop¬ 
ping  hydraulic  head  and  failure  time 


2'3 


Mgure  7. 


theoretical  explanation  restrict  the  applicability  of  those  empirical  rela¬ 
tions  and  indicate  the  need  for  a  more  thorough  and  detailed  analysis  of 
breaching  mechanisms. 


Mathematical  Modeling  of  Flood  Routing 


St.  Venant  system  of  equations 

35.  Propagation  of  a  flood  wave  through  the  receiving  channel  and  flood- 
plain  can  ba  successfully  described  by  means  of  unsteady,  incompressible, 
free-surface  hydrodynamic  equations.  More  specifically,  they  compose  the 
so-called  St.  Venant  system  of  equations  expressed  as 


3A  3Q 
3t  3x 


q^  (Continuity) 


(10) 


3Q  +  3_ 
3t  3x 


S  )  =  0  (Momentum) 
<V 


(ID 


where 

A  =  wetted  cross  section 

t  and  x  =  time  and  distance  coordinates,  respectively 
Q  =  water  discharge 

q  =  lateral  inflow 
o 

y  =  water  depth 

Sj.  =  energy  loss  gradient 

S  =  slope  of  the  channel 
o 

To  determine  the  energy  gradient,  either  Chezy's  or  Manning’s  friction  rela¬ 
tion  can  be  applied.  For  completeness  of  the  problem,  both  the  initial  and 
boundary  conditions  must  be  provided. 

36.  The  St.  Venant  system  of  equations  is  a  nonlinear  partial  differ¬ 
ential  system  of  the  hyperbolic  type  for  which  no  general  analytical  solution 
is  known.  Solution  of  that  system  can  be  obtained  only  by  means  of  three  main 
numerical  techniques:  the  characteristics,  the  finite  differences,  and  the 
finite  elements.  Each  method  is  described  in  the  following  paragraphs. 

37.  Characteristics  technique.  The  main  feature  of  this  method  is  the 
transformation  o f  the  original  partial  differential  system  of  two  equations 


into  an  ordinary  differential  system  of  four  equations.  This  is  possible 
because  the  system  is  hyperbolic.  The  characteristics  can  be  defined  as 
propagation  paths  of  a  geometric  or  physical  disturbance.  For  a  channel  of 
constant  width  and  zero  lateral  inflow,  Equations  10  and  11  can  be  written  as 
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dt 
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Combining  Equations  12  and  13  and  the  total  differentials  (du,dy)  yields 
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Equation  14  has  a  defined  solution  if  and  only  if  (Abbott  1966): 


u  4-  c 
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where 


C-  =  wave  characteristics 

J~  =  Riemann's  quasi-invariants 

c  =  wave  celerity 
o 

38.  Thus,  the  system  of  Equations  12  and  13  has  been  transformed  into  the 
system  of  Equations  15-18.  The  new  system  can  be  solved  graphically  (Schon- 
feld  1951) ,  semigraphically  (Chow  1959) »  or  numerically.  The  numerical 
solution  is  based  on  the  finite-difference  techniques.  The  solution  can  be 
obtained  either  on  a  characteristics  grid  (Figure  8)  in  explicit  form  (Faure 
and  Nahas  1961)  or  implicit  form  (Amein  1966)  ,  or  on  a  fixed  grid  (Figure  9) 
in  explicit  form  (Stoker  1957)  or  implicit  form  (Mozayeny  and  Song  1969)  . 

39.  Finite-difference  technique.  The  main  feature  of  finite-difference 
techniques  is  approximation  of  the  derivatives  in  the  governing  equations  by 
truncated  Taylor  Series  so  that  the  solution  is  obtained  on  nodal  points  of  a 
rectangular  x-t  fixed-grid  system.  The  solution  proceeds  from  time  step  j 

to  time  step  j+1  .  If  the  computation  advances  by  solving  a  single  equation, 
the  numerical  scheme  is  explicit.  If  the  computation  requires  the  solution  of 
a  system  of  equations,  the  scheme  is  implicit.  Explicit  schemes  were  sug¬ 
gested  by  Isaacson,  Stoker,  and  Troesch  (1958),  by  Courant,  Freidrichs,  and 
Lewy  (1967),  by  Lax  and  Wendroff  (I960,  1964),  and  by  Dronkers  (1964). 

Implicit  schemes  were  given  by  Preissman,  Vasilier,  and  Abbott  (Mahmood  and 
Yevievich  1975)  and  by  Dronkers  (1969). 

40.  Finite-element  technique.  In  this  method  the  solution  domain  is 

subdivided  into  a  number  of  subdomains,  the  finite  elements,  and  for  each  ele- 

(e) 

ment  the  unknowns  x  are  approximated  in  discrete  form  as 

/  •.  m 

X(S)  =  £NiX.  (19) 


where 

N.  =  shape  functions 

l 

x^  =  value  of  the  unknowns  on  the  nodal  points 
m  =  number  of  nodes  of  each  element 

Substitution  of  the  approximate  solutions  (Equation  19)  into  the  governing 
equations  produces  an  error  that  is  minimized  either  by  means  of  variational 
calculus  or  by  the  more  general  method  of  weighted  residuals  (Finlayson  1972) . 
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Figure  8.  Two-dimensional  characteristics  grid 
(L  =  left,  P  =  point  of  determination,  R  =  right) 


Figure  9.  Characteristics  on  a  rectangular  fixed 
grid  (A  and  B  =  arbitrary  grid  points,  M  =  mid¬ 
point,  j  =  discrete  time,  i  =  discrete  longitudi¬ 
nal  space) 


In  chat  way,  a  local  algebraic  equation  is  derived  for  each  element.  After 
assemblage  of  all  local  equations  into  a  global  system,  solution  is  obtained 
by  solving  the  system  and  determining  the  values  of  variables  at  each  nodal 
point  (Zienkiewicz  1971).  Depending  on  the  form  of  shaping  functions,  the 
numerical  method  can  be  either  a  hybrid  finite  element-finite  difference 
scheme  or  a  pure  finite-element  scheme.  More  specifically,  if  =  N^(x)  , 
then  the  finite  element  discretization  is  done  only  for  the  space  coordinates, 
while  the  solution  marches  in  time  by  a  finite-differences  algorithm  (Taylor 
and  Davis  1973).  If  =  N.(x,t)  ,  then  the  solution  is  based  entirely  on 
finite-element  technique  (Scarlatos  1982). 

Simplified  approaches 

41.  Depending  on  the  physical  conditions,  the  St.  Venant  system  of  equa¬ 
tions  can  be  reduced  to  a  simpler  form  by  neglecting  one  or  more  of  the  terms 
in  the  momentum  equation  (Equation  11  or  13).  A  visualization  of  various 
approximations  can  be  given  as  follows: 

I - r -  Gravity 


3u  3u 

3F+U  3^  + 


+  g(Sf 
! _ 


Kinematic 


Diffusive 


1 — — — —  ■  . . . — 1 - -  Steady  Dynamic 

- - — . . . . . . ^ -  Full  Dynamic 

42.  The  advantage  of  these  approximations  is  primarily  the  simplification 
of  computational  requirements.  However,  the  physical  problem  itself  dictates 
which  one  of  the  approximate  forms  is  more  appropriate.  It  has  been  proven 
that  the  kinematic  wave  model  is  a  very  useful  technique  for  flood  routing. 

An  extensive  treatment  of  kinematic  wave  modeling  was  given  by  Sherman  and 
Singh  (1978,  1982).  Another  approach  to  flood  routing  is  the  Muskingum 
method,  where  the  dynamic  equation  is  replaced  by  an  empirical  relation 
between  water  storage  and  inflow-outflow  discharges  (Singh  and  McCann  1980). 
Classification  and  comparison  of  models 

43.  The  St.  Venant  system  written  in  the  form  of  Equations  10  and  11 
neglects  the  effects  of  wind  stresses,  atmospheric  pressure  differences,  and 
the  Coriolis  Force.  Knowledge  of  initial  and  boundary  conditions  is  also 
required.  Experience  with  the  full  dynamic  model  has  shown  that  it  can  yield 
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results  of  sufficient  accuracy,  but  the  solution  is  sensitive  and  sometimes 
leads  to  computational  instabilities.  On  the  other  hand,  simplified  models 
show  a  more  stable  solution  behavior,  and  they  produce  some  kind  of  results 
under  all  circumstances.  In  many  cases,  however,  these  results  are  very 
inaccurate  and  of  no  practical  use.  Precise  delineation  of  conditions  under 
which  simplified  models  can  be  successfully  applied  has  not  yet  been  achieved 
The  problem  of  defining  the  best  model  is  very  complicated  due  to  the  large 
number  of  variables  involved.  Additional  confusion  is  introduced  by  the 
special  features  of  the  numerical  solution  technique  itself.  When  local  and 
convective  acceleration  is  negligible,  the  diffusive  model  can  be  applied. 
Furthermore,  if  pressure  variation  is  small  in  comparison  to  gravity  and  fric 
tion  effects,  the  kinematic  wave  approach  is  suggested. 
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PART  III:  GRADUAL  DAM-BREAK  EVOLUTION  AND  FLOOD  PREDICTION 


Dam-Break  Evolution 


44.  Simulation  of  the  total  earth-fill  dam-breach  erosion  process  is  a 
combination  of  hydrologic  elements,  hydrodynamics,  sediment  transport  mechan¬ 
ics,  and  geotechnical  aspects.  The  real-life  problem  is  unsteady,  nonhomo- 
geneous,  nonlinear,  and  three-dimensional,  which  is  not  theoretically  well 
understood.  Mathematical  modeling  of  the  phenomenon  requires  idealization  of 
the  real-life  situation  so  chat  the  leading  physical  processes  can  be 
described  by  a  set  of  governing  equations.  Assumptions  on  which  the  governing 
equations  are  based,  the  ability  to  determine  certain  parameters  involved,  and 
accuracy  of  the  solution  algorithm  control  the  validity  of  the  model.  For 
practicality ,  there  is  always  a  trade-off  involving  complexity,  accuracy,  and 
efficiency  of  the  model. 

45.  Earth-fill  dam-breach  erosion  is  understood  intuitively  as  a  two- 
phase  wa  er-soil  interacting  system.  Water  from  the  reservoir  flows  through 
the  breached  section  of  the  dam,  causing  enlargement  of  the  breach  either  by 
erosion  or  sloughing.  The  process  continues  until  the  reservoir  is  emptied  or 
the  dam  resists  further  erosion.  In  the  following  sections,  each  component 
and  process  of  the  reservoir-dam  system  will  be  presented.  Assumptions  and 
simplifications  will  be  discussed  and  explained  through  physical  reasoning. 
Reservoir  water  mass  balance 

46.  The  volume  of  water  stored  within  the  reservoir  V  is  a  function  of 
the  reservoir  geometry.  Theoretically,  this  volume  can  be  estimated  as 

H 

V  =  /A  (H)dH  (20) 

0  s 

where 

H  =  reservoir  water  level  measured  from  a  reference  datum 

A  =  surface  water  area  within  the  reservoir 
s 

Equation  20  assumes  a  horizontal  water  surface  within  the  reservoir,  neglect¬ 
ing  any  possible  surface  profile,  which  is  for  practical  purposes  correct 
und~r  equilibrium  conditions.  When  the  dam  is  breached,  water  from  the  equi¬ 
librium  stage  within  the  reservoir  starts  to  accelerate  and  converge  toward 
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the  breach,  while  at  the  same  time  there  is  a  continuous  depletion  of  the 
water  volume  V  .  This  phenomenon  is  essentially  dynamic  and  is  controlled  by 
both  the  mass  continuity  and  momentum  balance  equations.  Due  to  comparatively 
small  velocities  within  the  reservoir  and  the  locality  of  the  dynamic  effects, 
the  rate  of  water  volume  depletion  can  be  described  by  a  single  mass  continu¬ 
ity  equation' as 


d V 
dt 


(21) 


where 

I  =  inflow 
o 

=  breach  outflow  discharge 

Q  =  outflow  over  the  crest  of  the  dam 
o 

Qgp  =  outflow  through  the  spillway  and  powerhouse  outlet 
The  time  derivative  of  the  water  volume  can  be  written  as 


dV  =  dV  dH  =  ,  dll 
dt  dH  dt  s'1  }  dt 


(22) 


where  V  is  the  reservoir  water  storage  capacity.  Combining  Equation  21  and 
22  yields 


A  (H)  jyS  =  I  -  Q  -  Q  _  q  (23) 

s  dt  o  b  ^o  sp 

47.  Inflow  discharge  Iq  includes  all  water  sources  such  as  riverine 
water,  watershed  runoff,  direct  precipitation,  and  ground-water  flow  into  the 
reservoir.  The  combined  effect  is  given  in  the  form  of  a  hydrograph  through 
statistical  evaluation  of  existing  data.  The  more  extensive  and  accurate  the 
data  set,  the  more  reliable  the  inflow  hydrograph.  In  case  of  limited  data, 
an  inflow  hydrograph  should  be  assumed  that  corresponds  as  well  as  possible  to 
the  expected  conditions. 

48.  Another  specified  variable  is  the  outflow  Q  .  Indeed,  the 

sp 

spillway  capacity  is  given  as  a  function  of  the  water  elevation  H  ,  while  the 
powerhouse  discharge  is  also  a  predetermined  function  of  water  jlevation  and 
time.  Knowledge  of  both  of  these  quantities  is  essential  for  efficient 
operation  and  management  of  the  dam,  so  they  are  always  accurately  specified. 
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49.  Before  construction  of  a  dam,  the  upstream  valley  that  will  serve  as 

the  artificial  lake  is  mapped  in  detail  to  determine  the  storage  capacity  of 

the  reservoir.  Therefore,  the  relation  A  =  A  (H)  is  a  known  function.  m 

s  s 

most  cases,  however,  instead  of  the  Ag  =  Ag(H)  relation,  an  equivalent  rela¬ 
tion  of  V  =  V(H)  is  provided  so  that  the  Ag (H)  function  can  be  obtained 
directly  as  the  tangent  at  any  point  of  the  V-H  curve. 

50.  Referring  to  Equation  23,  it  is  obvious  that  the  only  unspecified 

quantities  are  the  outflow  discharges  through  the  breach  and  over  the  crest  of 

the  dam.  If  those  quantities  could  be  expressed  in  terms  of  only  the  water 

elevation  H  ,  then  Equation  23  would  be  an  ordinary  differential  equation 

that  can  be  solved  easily.  However,  as  it  will  be  shown  in  the  next  section, 

breach  outflow  Q,  contains  another  unknown  variable,  the  breach  bottom 
b 

elevation  Z  ,  so  that  Equation  23  cannot  be  solved  directly.  A  schematic 
presentation  of  the  geometric  and  physical  quantities  of  dam-break  problems  is 
given  as  Figure  10. 


Figure  10.  Geometric  end  physical  characteristics  of 
earth-dam  failure 
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Hydraulics  of  flow  through 
the  breach  and  ov<,r  *~he  crest 

51.  Flow  through  the  breach  and  over  the  crest  of  the  dam  resembles  flow 

over  a  broad-crested  weir.  Since  there  is  no  information  for  unsteady  broad- 
crested  weir  flow,  steady-state  expressions  for  the  flow  will  be  used  in  this 
study.  This  is  justified  by  the  fact  that,  in  the  vicinity  of  the  breach, 
local  accelerations  are  much  smaller  than  convective  accelerations  as  the  par¬ 
ticles  start  moving  from  rest  toward  the  breach.  Therefore,  quasi-steadv  con¬ 
ditions  will  describe  the  phenomenon  fairly  well,  and  both  outflows  through 
the  breach  and  over  the  crest  of  the  dam  will  be  taken  as 

Qb  =  [Cjb  +  C2  (H  -  Z)  tan  6 ] (H  -  Z)  '  (24) 

Q0  =  C*(Bd  -  B)(H  -  Zo)3/2  (25) 

where 

★  * 

and  C2  =  dimensional  coefficients 

o  -  bottom  width  of  the  breach 
Z  =  bottom  elevation  of  the  breach 
0  =  angle  between  vertical  and  the  breach  side 
Bp  =  top  width  of  the  dam  (crest  length) 

B  =  top  width  of  the  breach 

Equation  24  corresponds  to  a  trapezoidal-sh3ped  breach,  while  Equation  25  cor¬ 
responds  to  a  rectangular-shaped  weir.  For  b  =  0  ,  Equation  24  describes  a 
triangular  breach,  and  for  0  =  0  ,  a  rectangular  one. 

* 

52.  In  the  case  of  a  rectangular  weir,  the  theoretical  value  for  can 

be  easily  derived  from  critical  flow  conditions  over  the  crest  as 


=  1 . 7b (H  -  Z)3/2 


(26) 
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where 


A^  =  wetted  cross  section  of  the  breach 

y  =  critical  depth 
c  * 

Therefore,  in  the  metric  unit  system,  C  =  1.7.  The  theoretical  value  for 

*  1 

the  C0  coefficient  in  the  same  unit  system  is  1.35.  In  practice,  those 
values  should  be  reduced  due  to  correction  for  velocity  of  approach  (Brater 
1959)  . 

■k  k 

53.  Further  reduction  of  the  values  of  coefficients  and  C0  might  be 

necessary  when  tailwater  effects  are  present,  i.e.,  when  flow  is  submerged 
(Figure  11).  In  that  case,  these  coefficients  are  modified  from  the  equation 


C 


*m 

1,2 


(27) 


where 


C 


*m 


modified  C  coefficient 
I  ,  Z 


1,2 

■k  k  k 

Cl,2  =  C1  °r  C2 

y^  =  water  depth  at  the  tailwater  section 
Equation  27  is  an  empirical  relation  and  implies  that 
submergence  over  hydraulic  head  is  less  than  0.67,  the 
negligible . 


54.  The  water  depth 


v 

'  o 


is  computed  from  Chezv's 


if  the  ratio  of  depth  of 
tailwater  effects  are 

equation 


Wo'"2* 


(28) 


or  Manning's  equation 


1  2/3  1/2 

Q  =  -  R"  S  A 
n  h  o 


(29) 


where 

C,  =  Chezv's  coefficient  of  friction 
h 

R,  =  hvdraulic  radius  at  the  tailwater  cross  section 
h 

n  =  Manning's  coefficient  of  friction 
Equations  28  and  29  are  transcendental  equations  with  respect  to  v  and 
require  a  t ri al-and-er ror  procedure  for  their  solution. 


Figure  11.  Submerged  flow  conditions 

55.  Combining  equations  23,  24,  and  25,  the  reservoir  water  volume  deple¬ 
tion  equation  reads 

dF  *  *  3/2 

A  (H)  —  =  I  (t)  -  [C.b  +  C0(H  -  Z)  tan  0](H  -  Z) 
s  dt  o  12 

-  C* ( B  -  B)(H  -  Z  )3/2  -  Q  (H ,  t )  ( 301 

ID  o  sp 

Equation  30  is  a  nonlinear  ordinary  differential  equation  with  two  unknowns: 
water  elevation  H  and  breach  bottom  elevation  Z  .  Those  two  unknowns  are 
interdependent  through  the  processes  of  outflow  discharge  and  breach  erosion. 
For  completeness  of  the  solution,  an  equation  that  describes  dam-erosion 
characteristics  should  be  derived. 

Flow  through  breach  on 
the  downstream  face  of  dam 

56.  The  main  erosive  force  is  water  flowing  at  high  velocities  over  the 
downstream  face  of  the  dam.  Although  the  flow  is  unsteady,  it  can  be  approxi¬ 
mated  by  quasi-steady-state  conditions  by  the  same  reasoning  used  for  the  flow 
over  the  crest.  According  to  experimental  data  of  Pugh  and  Gray  (1984'  ,  the 
flow  over  the  whole  top  section  of  the  breach  can  be  assumed  as  being  critical 
(Figure  12).  Therefore,  the  water  flow  over  the  downstream  face  of  the  dam 
will  be  supercritical,  reaching  normal  flow  conditions  after  passing  through 
an  S2  profile  (Figure  10). 

57.  When  local  accelerations  are  neglected,  the  momentum  equation 
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(Equation  13)  can  be  written  as 


°  Critical  dep 


Embankment 


Figure  12.  Flow  over  the  crest  breach  section  (after  Pugh 

and  Gray  1984) 


d 

dx  \  2 

V2gAb 


Making  use  of  Chezy's  friction  equation  for  and  after  some  mathematica 

manipulations,  Equation  31  yields 


dx  l1  . 3  o  ,2,2, 


For  steep  slopes.  Equation  31  should  be  corrected  as 


S°  ChAbRh 


cos  a,  -  - r 

b  .  3 


where  a  is  the  angle  of  the  downstream  face  of  dam  with  the  horizontal . 
Integration  of  Equation  32  requires  an  iterative  technique.  Since  flow  is 


supercritical,  the  integration  starts  from  the  upstream  boundary,  i.e.,  the 
critical  depth. 

Erosion  processes  and  sediment  transport 

58.  After  development  of  an  initial  breach  on  the  dam,  the  hydrodynamic 
forces  continue  to  enlarge  the  breach  by  eroding  the  soil  material.  Mechanics 
of  sediment  transport  is  a  scientific  discipline  that  has  been  developed  in  a 
semiempirical  form  mostly  for  the  case  of  alluvial  rivers.  Because  of  a  lack 
of  information  on  sediment  erosion  under  extremely  dynamic  conditions,  such  as 
those  occurring  during  an  earth-fill  dam  failure,  sediment  discharge  will  be 
estimated  by  a  conventional  method,  the  Einstein-Brown  bed-load  formula  (Brown 
1950)  .  Although  this  method  has  been  successfully  applied  for  prediction  of 
sediment  transport  in  alluvial  streams,  its  application  to  dam-erosion  dynam¬ 
ics  requires  extrapolation  beyond  the  range  for  which  experimental  data  exist. 
The  Einstein-Brown  formula  was  chosen  since  it  has  been  more  widely  tested 
than  any  other  method  (Simons  and  Senturk  1976).  Besides,  this  method  does 
not  depend  on  a  threshold  value  of  shear  stress  for  initiation  of  erosion, 
which  cannot  be  determined  easily. 

59.  Einstein-Brown  bed-load  formula.  The  basic  idea  of  the  Einstein- 
Brown  theory  is  that  initiation  and  cessation  of  sediment  motion  depend  on  the 
probability  that  relates  instantaneous  hydrodynamic  lift  forces  to  the  sub¬ 
merged  weight  of  a  particle.  Their  final  results  are  presented  in  the  dimen¬ 
sionless  expression 


where 

=  sediment  transport  rate  function 
T  =  inverse  of  Shield's  dimensionless  shear  stress 
Explicitly,  the  quantities  <J>  and  E  are  given  as 


and 


(34) 


J_  _  _ T _ 

•  "  (Y2  -  Y)Ds 


(35) 
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where 


*bw 

Ys  = 
= 


bed-load  discharge,  weight  per  unit  width 

specific  weight  of  soil 

constant 


36v 


8D^-’ 


(36) 


Y  = 

D  = 
s 

T  = 


specific  weight  of  water 
representative  size  of  bed  sediment 
bed  shear  stress 

specific  weight  of  submerged  soil 
kinematic  viscosity  of  the  water 


Usually,  D  is  taken  as  the  median  size  D...  ,  while  bed  shear  stress  is 
s  50 

estimated  as 


T  =  YVf  =  Y  72 

Lh 


(37) 


60.  The  functional  relationship  of  Equation  33  was  determined  using 
experimental  data.  A  plot  of  the  results  is  given  as  Figure  13.  As  shown  in 
this  figure,  when  1/T  >  0.09  ,  Equation  33  becomes 


$  =  40  (  | 


(38) 


At  this  point,  it  should  be  mentioned  that  due  to  high  shear  stresses  exper¬ 
ienced  in  the  dam-erosion  problem,  the  value  of  l/Y  will  be  much  higher  than 
the  limiting  number  of  2  given  in  Figure  13.  Therefore,  in  that  case,  an 
extrapolation  will  be  necessary. 

61.  Breach  bottom  erosion  rate.  Once  the  bed-load  discharge  q^w  has 
been  estimated,  the  rate  of  erosion  of  the  bottom  of  the  breach  can  be 
directly  calculated.  Indeed,  scouring  AZ  of  the  breach  during  time  interval 
At  ran  be  given  as 
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AZ 


qbwAt 


Yg ( 1  “  pH 


(39) 


where  p  is  the  soil  porosity. 

62.  Since  bed-load  discharge  depends  on  hydrodynamic  conditions  and  those 
conditions  change  from  critical  to  supercritical  flow,  erosion  processes  must 
be  considered  separately  for  the  breach  at  the  crest  and  the  downstream  face 
of  dam. 

Geotechnical  considera- 

tions  of  breach  slope  stability 

63.  During  the  erosion  processes  of  an  earth-fill  dam,  the  situation 
arises  where  breach  slopes  become  unstable.  This  happens  when  the  hydrody¬ 
namic  forces  associated  with  seepage  are  greater  than  the  soil  friction  and 
cohesion.  The  problem  can  be  successfully  analy’ed  by  the  contour  method 
(Chugaev  1964),  in  which  the  shearing  surface  is  assumed,  for  simplicity,  as  a 
single  plane  passing  through  the  toe  of  the  slope.  A  schematic  representation 
of  the  problem  is  given  as  Figure  14.  The  initial  water  table  is  the  horizon¬ 
tal  line  3-4.  Due  to  breaching  and  depletion  of  the  reservoir  water,  the 
water  surface  is  drawn  down  to  line  2-5,  which  will  create  a  horizontal  seep¬ 
age  force  that  along  with  gravity  forces  might  cause  failure  of  the  slope. 

The  main  advantage  of  the  contour  method  is  that  it  requires  knowledge  of  the 
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Figure  14.  Characteristics  of  slope  instability 

head  distribution  only  along  the  boundaries  of  the  sliding  wedge  and  not 
throughout  the  entire  wedge. 

64.  In  the  contour  method,  the  total  seepage  force  acting  on  the  wedge  is 
obtained  directly  from  the  hydraulic  head  distribution.  Let  1-7-8  be  the 
sliding  wedge.  The  piezometric  line  of  the  upper  part  of  the  wedge  is  repre¬ 
sented  by  line  4-3-6-2-5,  while  the  piezometric  line  for  the  shearing  surface 
is  given  by  straight  line  4-6-5.  Projecting  the  hydraulic  heads  on  l'-l" 
axis,  the  horizontal  component  F  of  the  total  seepage  force  is  proportional 
to  the  area  of  triangle  1  *  — 2 *  —3 * • 

65.  For  estimation  of  the  weight  of  the  wedge,  the  nonuniform  presence  of 
water  within  the  sliding  wedge  should  be  considered.  Indeed,  the  total  weight 
of  the  wedge  can  be  estimated  by  calculating  the  weight  of  the  saturated  soil 
as  well  as  the  buoyancy  effects  as  follows.  Section  3-4-7-8  is  composed  of 
dry  soil  (y  ) .  In  section  2-6-5,  negative  pressure  is  assumed,  so  that  the 
specific  weight  is  that  of  pure  water  but  with  a  minus  sign  (-y)  .  Soil  is 
saturated  in  section  3-6-4,  so  the  specific  weight  is 

Y  =  Y  +  PY  (40) 

1  s 

where  y  is  the  specific  weight  of  saturated  soil.  Finally,  the  soil  sec¬ 
tion  1-4-6  is  submerged,  with  specific  weight  y?  given  as 


40 


Y  =  Y  -  (1  -  p) Y 
l  s 


(41) 


The  total  weight  of  sliding  wedge,  G  ,  is  the  sum  of  the  four  separate  parts, 
and  for  a  wedge  of  unit  width  it  yields 

G  =  YsA(3-4-7-8)  ~  yA(2-6-5)  +  YlA(3-6-4)  +  Y2A(l-4-6)  (42) 

where  is  the  area  of  each  individual  section. 

66.  Stability  or  failure  of  the  breach  sides  depends  upon  the  balance  of 
forces  acting  on  the  wedge.  Those  forces  are  the  weight  of  the  wedge,  the 
seepage  forces,  the  internal  friction,  and  the  cohesion.  At  the  stage  of 
equilibrium,  the  force  balance  equation  yields  (Chugaev  1964) 

F  +  G  tan  (4  -  <f>)  =  Cx  [ 1  +  tan  £  tan  (t;  -  4>)  ]  (43) 

H  p 


where 

4  =  angle  between  the  shearing  plane  and  the  horizontal 
C  =  cohesion 

Xp  =  horizontal  projection  of  the  shearing  plane 
Failure  of  the  slope  occurs  when  the  right-hand  side  of  Equation  43  is  greater 
than  the  left-hand  side. 

Flood  Routing  by  the  Muskingum-Cunge  Method 

67.  Once  the  outflow  hydrograph  from  the  breach  is  known,  the  flash  flood 
can  be  routed  through  the  downstream  receiving  channel.  One  well-established 
technique  for  flood  routing  is  the  Muskingum-Cunge  method  (Ponce  and  Yevjevich 
1978).  This  method  is  based  on  a  linear  relation  between  inflow  I  ,  outflow 
0  ,  and  reach  storage  S  ,  given  the  form 

S  =  K[al  +  (1  -  a)0]  (44) 

where 

K  =  dimensional  coefficient 
a  =  weighting  factor 

Equation  44  is  coupled  with  the  volume  continuity  equation  written  as 
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(45) 


In  Equations  44  and  45,  the  function  I  is  known  either  from  the  flash-flood 
hydrograph  or  from  the  computations  of  the  adjacent  upstream  reach. 

68.  In  contrast  with  the  original  Muskingum  method  where  both  K  and  a 
are  constant  parameters,  in  the  Muskingum-Cunge  method,  K  and  a  vary 
according  to  the  expressions 


K  = 


Ax 

c 

o 


(46) 


and 


where 


(47) 


Ax  =  length  of  a  channel  reach 

c  =  wave  celerity 
o  J 

q  =  discharge  of  unit  width 

It  has  been  proven  that  application  of  this  routing  technique  can  give  results 
comparable  in  accuracy  to  the  application  of  the  diffusive  model  (Ponce  and 
Yevjevich  1978). 


Numerical  Solutions  of  the  Governing  Equations 


69.  Once  the  governing  equations  have  been  defined,  the  next  step  is  to 
determine  their  solution  algorithm.  Unfortunately,  most  of  the  equations  can¬ 
not  be  solved  analytically,  so  a  numerical  solution  is  required.  In  this  sec¬ 
tion,  emphasis  will  be  restricted  to  certain  independent  solution  techniques 
and  not  the  overall  dam-break  problem. 

Solution  of  the  water-profile  equation 

70.  For  the  solution  of  the  water-profile  relation  (Equation  32)  ,  the 
numerical  technique  suggested  by  Prasad  (1970)  will  be  used.  Let  the  flow 
profile  be  described  by  y  =  f (x)  .  Applying  the  trapezoidal  rule  of 
integration , 


dv 

dx 


i+1 


*i  + 


i+1 


+  ^ 

dx 


Ax 


(48) 
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where  the  subscript  i  refers  to  the  distance  along  the  channel  and  it 
increases  downstream. 

71.  Based  on  Equations  32  and  48,  the  two  unknowns  y  and  dy/dx  can  be 
computed  as  follows: 

Estimate  (dy/dx) |  from  Equation  32,  either  from  initial 
data  or  previous  calculation. 

Set  (dy/dx) | =  (dy/dx) I  as  a  first  approximation. 
Obtain  an  approximate  value  for  y^+^  from  Equation  48. 
Compute  a  new  value  for  (dy/dx) from  Equation  32  using 
the  obtained  in  step  3. 

If  the  new  value  of  (dy/dx) | is  not  very  close  to  the 
value  previously  assumed  or  computed,  then  repeat  steps  3-5. 
Otherwise,  proceed  to  the  next  integration  step  and  repeat 
the  whole  procedure. 

The  method  is  fast  and  accurate  and  can  be  programmed  very  easily. 

Solution  of  the  Muskingum-Cunge  equation 

72.  Combining  Equations  44-47  and  setting  them  in  finite-difference  form, 
after  some  manipulations ,  results  in  the  following  equation: 

0j+1  =  C1Ij  +  C2Ij+1  +  C30j  (49) 

where  the  upper  index  j  refers  to  the  time  step  and  C^  ,  ,  C^  are 

numerical  coefficients.  The  space-time  discretization  of  the  Muskingum-Cunge 
method  is  shown  in  Figure  15.  From  this  figure  it  is  evident  that  the  outflow 
of  a  specific  section  is  inflow  for  the  downstream  adjacent  section. 

73.  The  coefficients  C^  ,  ,  and  C^  can  be  evaluated,  respectively, 

from  the  following  relations: 

1  +  C  -  C 

C  =  --  4  5 

1  1  +  c.  +  c, 

4  5 

-1  +  c4  +  c5 

C2  1  +  C  +  C 
4  5 

1  -  c  +  c 

c  - - - 

3  1  +  C.  +  C, 

4  5 


(50) 

(51) 

(52) 


Step  1. 

Step  2. 
Step  3. 
Step  4. 

Step  5. 
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Figure  15.  Space-time  discretization  for  the 
Muskingum-Cunge  method 


in  which  C.  and  Cc  are  defined  as 
4  5 


C  =  c 
4  o\  Ax 


C,  =  — 7T 
5  c  Ax 
o 


The  time  step  At  is  usually  taken  as  constant.  Both  and  C,_  have 

physical  significance,  being  a  ratio  of  celerities  and  dif fusivi ties , 
respectively. 

74.  For  the  estimation  of  these  coefficients,  it  is  necessary  to  deter¬ 
mine  the  wave  celerity  c  and  the  unit  width  discharge  q  for  each  computa 

o 

tional  cell.  The  values  of  c  and  q  are  defined  as 


where  T  is  the  top  width  of  the  channel  netted  cross  section.  To  compute 
coefficients  C,  and  Cc  ,  both  c  and  q  are  obtained  directly  as  a 

HD  O 

three-point  average  of  their  values  at  points  (i,j),  (i,j+l),  and  (i-*-l.j). 
This  method  has  been  proven  sufficiently  accurate  in  the  simulation  of  flood 
flows  (Ponce  and  Yevjevich  1978). 

75.  In  many  cases,  especially  when  dealing  with  trapezoidal  cross 
sections,  the  situation  arises  when  the  roots  of  an  implicit  algebraic  func¬ 
tion  y  =  f(x)  should  be  determined.  The  most  commonly  used  technique  for 
that  purpose  is  the  Newton-Raphson  iteration  algorithm,  given  as  follows: 


f(x.) 

Xi+1  Xi  f ' (x.) 

l 


where  i  is  the  iteration  index  and  f'(x)  is  the  first  derivative.  The 
method  is  very  efficient  and  converges  rapidly. 

Fixed-point  iteration  algorithm 

76.  In  certain  cases,  it  is  very  convenient  to  use  a  more  simplistic 
iteration  algorithm  such  as  the  fixed-point  scheme  instead  of  the  Newton- 
Raphson  technique.  A  graphical  description  of  that  scheme  is  presented  as 
Figure  16. 


Figure  16.  Fixed-point  iteration  algorithm 


PART  IV:  COMPUTER  MODEL  FOR  BREACH  EROSION  OF  EARTH-FILL  DAMS 


77.  The  Breach  Erosion  of  Earth-Fill  Dams  (BEED)  computer  model  is  a 
mathematical  model  developed  for  predicting  the  hydrograph  of  a  flash  flood 
due  to  gradual  dam  failure.  The  structure  of  the  model  is  based  on  the  quan¬ 
titative  and  aualitative  physical  nrinr’ples  described  in  Parts  IT  and  TTT. 

The  solution  procedures  and  algorithms  of  the  model  are  relatively  simple  and 
can  be  used  in  both  microcomputers  and  mainframe  computer  systems. 

Physical  Description  of  BEED 

78.  Before  presenting  a  quantitative  description  nj  che  BEED  model,  it  is 
important  to  examine  the  conceptual  framework  of  the  model  and  to  discuss  its 
physical  reasoning  and  consistency  as  well  as  its  applicability  and 
limitations . 

79.  The  model  will  be  Ge>’ for  a  homogeneous  dam  with  different  but 
uniform  slopes  for  the  upstream  and  downstream  faces.  Physical  and  geometric 
cnaracteristics  of  the  dam  and  its  surroundings  should  be  specified.  The 
model  neglects  the  triggering  mechanism  of  failure  and  can  simulate  the  phe¬ 
nomenon  only  when  a  small  breach  has  been  developed  at  the  crest  of  the  dam. 
The  size,  shape,  and  location  of  this  initial  breach  should  be  provided  as 
initial  conditions.  Unfortunately,  the  selection  of  such  conditions  is  based 
entirely  on  engineering  judgment  and  not  on  quantitative  information.  For 
convenience,  a  rectangular  initial  breach  shape  with  specified  depth-over- 
width  ratio  can  be  assumed. 

80.  Once  the  initial  breach  has  occurred,  water  from  the  reservoir  starts 
flowing  through  the  breach,  causing  enlargement  of  the  breach  and  erosion  of 
the  downstream  face  of  the  dam.  The  erosion  is  restricted  to  a  channel  of  the 
same  top  width  as  the  breach  at  the  crest  of  the  dam.  However,  the  erosion 
processes  occurring  on  the  crest  and  on  the  downstream  face  are  considered 
separately  because  the  water  velocities  are  much  higher  down  the  face  of  the 
dam  than  they  are  over  the  crest;  consequently,  the  downstream  face  erodes 
much  faster. 

81.  The  enlargement  of  the  top  breach  follows  a  similar  pattern  except  in 
the  case  where  sloughing  due  to  slope  instability  occurs.  At  this  point,  it 
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should  be  emphasized  that  the  model  incorporates  only  rectangular-  or 
trapezoidal-shaped  breaches. 

82.  Sloughing  effects  are  considered  only  for  the  crest  breach;  the  shape 
of  the  eroded  channel  on  the  downstream  face  is  adjusted  within  the  model  to 
that  of  the  crest.  If  the  normal  flow  conditions  at  the  downstream  face 
acquire  ccnciderable  ’istance  to  be  developed,  the  slope  should  be  subdivided 
into  a  number  of  reaches  so  that  the  erosion  in  each  can  be  estimated  sepa¬ 
rately.  This  will  result  in  a  steeper  downstream  slope.  In  most  cases,  how¬ 
ever,  normal  flow  conditions  are  established  in  a  very  short  distance  so  that 
the  effect  of  the  S2  water  prof i Ip  can  be  neglected  and  the  erosion  of  the 
downstream  face  can  be  assumed  ao  being  uniform. 

83.  When  the  flow  conditions  are  submerged  (i.e.,  when  tailwater  effects 
are  present) ,  the  assumption  is  made  that  erosion  occurs  only  on  the  top 
breach  and  not  on  the  downstream  face  while  breach  outflow  discharge  is 
reduced. 

84.  Another  characteristic  feature  of  the  program  is  that  when  the 

upstream  and  downstream  slopes  of  the  dam  meet  at  a  single  point  S  ,  a  sudden 
predetermined  mass  erosion  is  considered  so  that  a  new  top  horizontal  breach 
channel  of  length  is  established.  The  upstream  face  slope  of  the  dam 

remains  unaltered  during  the  failure  processes.  The  program  continues  the 
simulation  until  either  the  reservoir  is  emptied  of  water  or  the  d— n  resists 
any  further  erosion. 

! 

I 

Solution  Algorithm  of  BEEP 

j  85.  The  BEED  model  is  designed  to  estimate  dam-breach  erosion  processes 

to  predict  outflow  discharge  and  to  route  it  through  the  receiving  channel. 
Because  of  the  implicit  form  of  the  governing  equations,  the  solution  algo¬ 
rithm  is  iterative.  However,  practical  experience  indicates  that  convergence 
is  achieved  after  few  iterations. 

I  86.  The  first  step  in  the  BEED  model  is  definition  of  the  geometric  and 

physical  features  of  the  system.  These  are  given  separately  for  the  dam,  the 
reservoir,  and  the  downstream  channel.  At  the  same  time,  all  preliminary  com¬ 
putations  are  executed,  while  initial  conditions  are  specified. 

87.  More  specifically,  dam  dimensions  are  provided  along  with  soil  char¬ 
acteristics  such  as  specific  weight,  particle  diameter  size,  cohesion, 
i 
i 


47 


internal  friction  angle,  and  surface  roughness.  Functional  relationships  for 
the  spillway,  powerhouse  outlets,  and  the  reservoir  capacity  are  also  speci¬ 
fied.  Description  of  the  downstream  receiving  channel  is  given  through  defi¬ 
nition  of  shape,  size,  and  roughness  of  a  certain  number  of  reaches,  which  may 
be  subdivided  into  smaller  segments  by  linear  interpolation  techniques. 
Finally,  tne  size  and  shape  of  the  initial  breach  as  well  as  the  initial 
hydraulic  head  are  specified  a  priori.  Having  all  the  required  information 
available,  the  model  proceeds  by  estimating  reservoir  water  level,  breach  bot¬ 
tom  elevation,  and  outflow  discharge,  which  is  subsequently  routed  downstream. 
88.  The  main  variables  of  the  problem  are  reservoir  water  level  H  and 

breach  bottom  elevation  Z  .  However,  the  breach  outflow  discharge  Q,  is 

b 

used  as  an  additional  variable  during  the  iteration  processes.  For  the  solu¬ 
tion,  Equation  23  is  discretized  as 


/ H . , .  -  H.\ 

A  (  — -  - -  =  4  (I  +  I  -  Q.  -  Q.  -  Q 

s . , . \  At  12  o .  ,,  o.  b.,.  b .  o.,. 

1+1X  /  3+1  1  1+1  1  1+1 


-  Q  -  Q  -  Q  ) 
o.  sp. , .  sp. 

1  1+1  rl 


(58) 


and  then  written  as 


„  +  C.  -  V,  -  V  -  V,  -  V  -  V,,  -  V 

H  =  H  .  +  —A — — - i - Hi - _1 - 111 - i - fill - lii  (59) 

l+l  l  ZA 

Si+1 


where  i  is  reterred  to  known  time  t  and  i+1  .s  referred  to  new  time 
t+At  .  In  Equation  59,  the  quantities  H  and  Z  are  involved  implicitlv. 

89.  The  computational  steps  for  estimation  of  the  variables  are  as 
follows : 

Step  1.  Set  H.,,  and  Z.  ,  equal  to  values  obtained  from  the 
l+l  l+l 

previous  time  step  nr  initial  conditions  (  denotes  uncor¬ 
rected  value)  . 

Step  2.  Compute  outflow  discharge  from  Equations  24  and  25 

(  denotes  estimated  value). 

Step  3.  Check  for  tailwnter  effects  and,  if  needed,  correct  > 

according  to  Equation  27. 
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Step  4.  Compute  sediment  transport  and  rate  of  erosion  from 

Equations  33,  34,  35,  and  39.  Then  estimate  ( ~~ 

denotes  corrected  value) . 


Step  5.  If  [  ^  i-1-1  ~  ^i+~'  ^  very  small,  proceed  to  the  next  step. 

Otherwise,  set  Z.  .  =  Z.  ,  and  return  to  step  2. 
l+l  l+l 

Step  6.  Compute  H  ^  from  Equation  59. 

Step  7.  If  iH.  ,  -  H.  ,  is  not  very  small,  proceed  to  the  next 
l+l  l+l 

step.  If  it  is  the  first  iteration,  go  to  step  9  or  return 
to  step  2. 

Step  8.  Set  =  ^q+q  ’  c^iec'<;  for  tailwater  effects,  and  return 

to  step  6. 

Step  9.  Check  for  slope  stability  using  Equation  43. 

Step  10.  Adjust  dam-breach  dimensions. 

Step  11.  Compute  total  outflow  discharge. 

Step  12.  If  hydraulic  head  h  =  H  -  Z  is  zero,  proceed  to  the  next 
step;  otherwise,  return  to  step  1. 

Step  13.  Estimate  the  Muskingum-Cunge  coefficients  from  Equations 
50-54. 


Step  14.  Route  the  flood  according  to  Equation  49. 

90.  The  solution  algorithm  of  the  BEED  model  is  represented  in  flowchart 
form  in  Figure  17.  The  effects  of  nonuniform  flow  over  the  downstream  face  of 
the  dam  were  neglected  since  normal  flow  conditions  were  attained  in  a  very 
short  distance  due  to  high  slopes. 
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Figure  17. 


Flowchart  of  BKFD  computer  model  (Continued) 


Figure  17.  (Concluded) 


PART  V:  ANALYTICAL  SOLUTIONS  OF  GRADUAL  DAM  EROSION 

91.  A  better  insight  of  the  physical  processes  and  the  significance  of 
the  controlling  parameters  of  the  gradual  failure  of  a  dam  can  be  obtained 
through  analytical  expressions.  Unfortunately,  the  governing  equations  are 
very  complicated,  so  that  a  general  analytical  approach  is  not  possible.  The 
degree  of  complexity  of  the  system  can  be  drastically  reduced  by  making  proper 
simplifying  assumptions  and  by  lumping  a  number  of  physical  parameters  into  a 
form  of  constant  coefficients.  In  that  case,  closed-form  solutions  are  feasi¬ 
ble.  The  analysis  in  this  section  is  based  primarily  on  the  principles  pre¬ 
sented  in  Part  III. 

92.  Assuming  that  the  inflow  into  the  reservoir  is  of  a  much  smaller 
order  of  magnitude  than  the  breach  outflow  discharge  and  neglecting  the  spill¬ 
way  and  powerhouse  outflow,  the  mass  continuity  relation  (Equation  23)  becomes 

A  ( H )  77  =  -Q.  (60) 

s  d  t  b 

Furthermore,  if  A^  is  independent  of  H  (i.e.,  prismatic  reservoir)  and  the 
outflow  is  taken  as 


Qb  =  uAb 


(61) 


then  Equation  60  yields 


s  dt 


(62) 


where  A^  can  be  rectangular,  triangular,  or  trapezoidal.  From  Equations  24 
and  25  it  is  evident  that  the  water  velocity  at  the  breach  can  be  estimated  as 


u  =  ax(H  -  Z) 


(63) 


where  and  are  proper  coefficients.  Combination  of  Fquations  62  and 

63  gives  a  single  equation  with  two  unknown  quantities:  the  water  depth  H 
within  the  reservoir  and  the  elevation  of  the  bottom  of  the  breach  Z  . 
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Therefore,  an  additional  equation  is  required  for  the  solution  of  the  problem 
that  should  be  obtained  from  the  mechanics  of  sediment  erosion. 

93.  It  is  known  that  the  rate  of  erosion  is  a  function  of  the  bottom 
shear  stresses  or,  equivalently,  the  square  of  the  mean  water  velocity. 
Mathematically,  this  can  be  expressed  as 


dZ  B2 

dt  =  "V 


(64) 


where  and  are  proper  coefficients.  Depending  on  the  value  of  the 
exponents  and  62  »  the  system  of  Equations  62  and  64  can  be  linear  or 
nonlinear.  For  completeness  of  the  problem,  the  initial  conditions 


H  =  H  and  Z  =  Z  at  t  =  t 
o  00 


(65) 


should  be  provided. 

Rectangular  Breach 

94.  For  rectangular  breach  of  constant  width  b  ,  the  cross  section 
is  given  as 


=  b(H  -  Z) 


(66) 


This  implies  that  the  breach  enlarges  only  in  the  vertical  direction. 

Linear  case 

95.  If  the  rate  of  erosion  is  a  linear  function  of  the  velocity  (6  =  1) 1 

then  the  whole  problem  is  linear.  Combining  Equations  62  and  66  yields 


A  ^  =  -ub(H  -  Z) 
s  d  t 


(67) 


Dividing  Equation  67  by  Equation  62  and  rearranging, 


dH 

dZ 


b 


a?A^ 


(H  -  Z) 


(68) 
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and  setting  h  =  H  -  Z  ,  then  dH/dZ  =  dh/dZ  +  1  .  Therefore, 


dH  _  b 

dZ  a~A 
2  s 


h  -  1 


(69) 


Separating  the  variables  and  integrating  gives 

ln  (z^r  -  9  ■ z  +  ci  <™> 

where  is  an  integration  constant.  From  the  initial  conditions  (Equa¬ 

tion  65)  ,  Cj  is  estimated  as 


a  A 

r  2  s  i 

ci =  —  ln 


b (H  -  Z  ) 

0  o  _ 

a0A 
2  s 


-  Z 


(71) 


Substituting  Equation  71  into  Equation  70  gives  (after  some  manipulations) 


a„A  /  a  A 

H  =  Z+—  +\Ho  '  Zo  6XP 


-V  (7.  -  Z) 

a„A  o 


'  s 


(72) 


Equation  72  prescribes  the  elevation  of  t:he  breach  bottom  Z  as  a  function  of 
the  water  height  H  and  breach  characteristics. 

96.  It  is  desirable,  however,  to  have  Z  as  a  direct  function  of  time. 
Specifying  coefficients  and  8^  as  =  V®  and  8^  =  1/2  ,  Equation  64 
becomes 


dZ 

d 


\  =  -a2>jg(H  -  Z) 


(7  3) 


Combining  Equations  72  and  73  yields 


dz  )  U"A 


dt 


=  -v  j-H +  \H 


o 


a»A  \ 

7  2  s 

Z  -  — r—  exp 
o  b 


(z  -  Z) 

cx„A  c 


1/2 


(74) 


By  separating  the  variables,  Equation  74  is  written  as 


Equation  80  specifies  the  progression  of  breaching  in  time. 

97.  Commenting  on  the  assumptions  made  for  the  derivation  of  Equation  80, 

the  one  that  deviates  mostly  from  reality  is  =  1  .  From  experimental 

results,  that  exponent  ranges  from  4  to  6  (Laursen  1956).  Numerically,  this 
can  be  partially  corrected  by  adjustment  of  the  coefficient  .  Another 

critical  point  is  the  assumption  of  constant  width  b  ,  which  is  very  unlikely 
to  occur  in  the  case  of  a  ±arge  dam. 

Nonlinear  case 

98.  In  this  case,  ^  ^  1  >  i-e.,  erosion  is  a  nonlinear  function  of  the 
velocity.  Dividing  Equation  62  by  Equation  64, 


dH  \ 

dZ  “  VsV  . 


(81) 


Combining  Equation  81  with  Equations  63  and  66, 


dH 

dZ 


a.A 
2  s 


1-8  1+6. (1-6, ) 

al  (H  -  Z) 


(82) 


Setting 


1-6, 


ba. 


A3  = 


a„A 
2  s 


(83) 


and 


A4  =  1  +  61(1  -  62) 


(84) 


Equation  82  reduces  to 


dH  7 \  4 

dZ  =  A3(H  *  Z) 


Equation  85  can  be  transformed  to 


-(1/A  )  A 

A3  fz+1=W 


(85) 


(86) 
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i/a4 

where  W  =  (H  -  Z)A^  .  Separation  of  the  variables  and  integration  of 
Equation  86  yields 

1/A 

=  -A3  +  CT  (87) 

The  left-hand  side  of  Equation  87  is  the  Bakhmeteff  function. 

99.  A  closed-form  solution  is  feasible  if  the  coefficient  A^  is  prop¬ 
erly  defined.  Assuming  83  =  2  ,  =  yjg  ,  and  6^  =  1/2  ,  then 


A 


1 


a2Asg 


1/2 


(88) 


and 


A  =i 
2  2 


(89) 


so  that  Equation  87  becomes 


f  dW 
J  1  -  w1 


72=-A3Z  +CI 


(90) 


By  setting  =  W  ,  Equation  90  can  be  easily  integrated  to  give 


W  =  In  [ 1 


-  \)  -  (~r) 


Z  +  C, 


(91) 


or  by  inserting  the  transformations  back. 


8  ,/u  ■?  _i_  i  _  1 1  8  / 


a2Asg 


172 


v'H  -  Z  +  In  1  - 


a2Asg 

,2 


172  ’ 


H  -  Z 


Ma2As)Z  +  Cl 


(92) 
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From  the  initial  conditions  (Equation  65)  ,  the  integration  constant  C^.  can 
be  estimated  as 


ci  ■ 


«2As8 


mf^o ♦ 1"  (‘  - 


Z_o  /_b_\2 

2g  U2AS 


Substitution  of  the  integration  constant  into  Equation  92  gives 


(93) 


a2V 


1/2 


/H  -  Z  -  /H  ~ 


Z  |+  In 
o 


.  1/2  7x 1 /2 

“2Agg  (H  -  Z) 

b  “  77V2 

m  (Ho  v 


Vs8 


(94) 


M=k  1 


(Z  -  zo) 


Equation  94  describes  breach  erosion  in  terms  of  the  hydraulic  head  H  -  Z  . 

100.  To  establish  the  variable  Z  as  an  explicit  function  of  time,  an 
expression  for  the  quantity  H  -  Z  should  be  derived.  For  that  purpose. 
Equation  64  is  subtracted  from  Equation  62  while  both  Equations  63  and  66  are 
used , 


HH-)=.!!i(H-z)1+b  +  v^H.z)V2 

s 


Specifying  the  variables  ,  8^  ,  and  B7  as  before, 


d(H  -  Z)  b  r  „.3/2  ,  7. 

=  ~  T~\g  (H  -  Z)  +  a  g(H  -  Z) 


dt 


Separating  the  variables  and  setting 


(95) 


(96) 
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w 


2 


_ _ b_ 

a2A2g 


Z 


then  Equation  96  becomes 


dW, 


gao 


or 


W2(l  -  W2)  2 


^  -  r>2  ’ 


dt 


(97) 


g“i 


dt 


(98) 


Integration  of  Equation  98,  determination  of  the  integration  constant,  and 
substitution  of  the  original  variables  provides 


a.A  /g(H  -  Z  ) 
2  s  o  o 


H  -  Z 


v^H  -  Z  -  (  b/H  -  Z”  -  a. A  v'g 


2** 

s 


g  exp  l- 


“2g 


(99) 


Having  the  e:  <-^sion  for  the  hydraulic  head  (Equation  99),  the  progression  of 
the  breaching  can  be  directly  estimated  from  Equation  94,  The  nonlinear  case 
is  an  improvement  of  the  linear  one  because  it  approximates  the  rate  of 
erosion  with  a  quadratic  velocity  function. 

Triangular  Breach 

101.  For  triangular  breach  geometry,  the  cross-sectional  area  A^  is 
given  as 


Ab  =  s(H  -  Z)2  (100) 

where  s  is  the  side  slope  (1V:sH).  The  assumption  of  constant  side  slope 
implies  that  the  breach  will  enlarge  in  a  similarity  pattern. 
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Linear  case 


102.  Again,  fcr  the  linear  case,  82  =  1  -  Combining  Equations  62  and  100 
and  dividing  by  Equation  64  yields 


=  — f-  (H  -  Z)2 
dZ  a2A2 


(101) 


Equation  101  can  be  rewritten  as 


d<H  -  Z)  +  1  -  -4-  (H-Z)2 


dZ 


a„A 
2  s 


(102) 


Setting  n  =  H  -  Z  in  Equation  100,  separating  the  variables,  and  using 
partial  fractions  gives 


dh 


dh 


M  2 


=  -2  dZ 


,a2As 


h  1  + 


a  A 
l  2  s 


(103) 


After  integration,  Equation  101  read« 


In 


1  +  ( — 
\a2As 


1/2 


a/2 


1  - 


,  ot  A 
\  2  s 


/  \l/2 

=  -2  4;)  z  +  c< 


(104) 


Estimating  the  integration  constant  C^.  from  the  initial  conditions  and  sub¬ 
stituting  back  to  Equation  104  gives  (after  some  algebraic  manipulations), 


103.  Equation  105  describes  the  changes  of  hydraulic  head  in  terms  of 
breach  bottom  elevation  Z  .  Therefore,  the  function  Z  =  Z(t)  must  be 
determined.  Equations  62,  63,  and  64  can  be  combined  to  give 

dH  dZ  s  1/2  2 . 5  ,  1/2,, ,  ,,1/2 

dt  "  dt  =  "  r  g  (H  “  2)  +  a28  (H  "  z) 

s 

Setting  h  =  H  -  Z  and  separating  the  variables  in  Equation  106  gives 


dh 


1/2/  s 


(=k 1,2 ' l) 


=  -«2g 


1/2 


dt 


Referring  to  page  72  of  Gradshteyn  and  Ryzhik  (1983),  the  integral  of 
Equation  107  is 


*1/4  ,1/2  . 

A.  -  h  1/2 

i  0  *.  -1  h 

In  —r~r, - r-rz  -  2  tan 


.1/4  .  ,1/2 
A,.  +  h 


1/4 


2  A^/4 

Vs  5 


-a,g1/2t  +  C, 


where  A,.  is  given  as 


A5  = 


ct-A 
2  s 


By  determining  the  integration  constant  C  ,  Equation  108  yields 


2  tan 


-1  (H  -  Z) 


1/2 


o„A 
2  s 


1/4 


-2 


3/4  1/4 
s 


+  In 


a„A 

2  s 


1/4 


a„A 
2  s 


1/4 


1/2  1/2 

-  (H  -  Z  )  '  .  (H  -  Z  ) 

o  o'  „  _  -1  o  o 

-  -  2  tan  - 


+  (H  +  Z  ) 
o  o 


1/2 


'jt„A 


1/2 


(106) 


(107) 


(108) 


(109) 


(110) 
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Equations  105  and  109  can  be  combined  to  determine  the  erosion  rate  of  breach 
bottom  or  the  depletion  of  reservoir  water. 

Nonlinear  case 

104.  For  the  nonlinear  case,  the  erosion  exponent  will  be  taken  as 
$2  =  2  ,  and  the  discharge  exponent  as  6^  =  1/2  .  Then,  dividing  Equation  62 
by  Equation-  64, 


dH 

dZ 


a2Asg 


1/2 


(H  -  Z) 


1.5 


(111) 


Setting  h  =  H  -  Z  and  separating  the  variables  in  Equation  ill  gives 


dh 


■i  ♦  h3/2 


=  dZ 


Vs6 


1/2 


Equation  112  can  be  easily  transformed  to 


A?/3 

W  dW  % 


1  -  w 


dZ 


where 


(112) 


(113) 


W 


.1/3,  1/2 
=  ^  h 


and 


a2Asg 


1/2 


(114) 


Referring  to  page  64  of  Gradshteyn  and  Ryxhik  (1983),  Equation  113  is  inte¬ 
grated  as 


I ln 


(i  -  wr 

1  +  w  +  w2 


1  -1  /  2W  +  1 

tan 


,1/2 


,1/2 


2/3 


2  z  +  cx 


or 


I  ln 


M73*1'2)2  ’ 

1  -1 

1/31/2 

2A^ +  1 

"  7  *1/3.  1/2  7  2/3. 

1  +  A,  h  +  A^  h 

I/O  Can 

3 

3 1  /  2 

2/3 


Z  +  C, 


Determining  the  integration  constant  from  initial  conditions,  the 

hydraulic  head  is  specified  as  a  function  of  Z  as 


ln 


!  +  Aj/3h1/2  +  Aj'V 
(1  -  A^/3h1/2)2 


z( 31 /2)  tan  1 


\ 


2A3/3h1/2  +  1 

TT72 


2/3 

3A^  (Z  -  Z)  +  ln 


1/31/2  2/3 

1  +  a£  V'  +  A^/3h( 


(i  -  Ay3hi/v 

D 


/  I  /  a  1  /  a 

'  V'  *•  +  1 


105.  To  determine  the  variable  Z  explicitly,  another  equation  is 
required.  For  this  purpose,  Equation  62  is  subtracted  from  Equation  61 


dH  _  dZ 
dt  dt 


dH  dZ  _  sg 


1/2 


(H  -  Z)2*3  I  a.,g(H  -  Z) 


By  separating  the  variables,  Equation  118  can  b«  written  as 


(115) 


(116) 


(11?) 


(1181 


6  3 


dh 


(119) 


,  ,,  A  ,3/2. 
h  ( 1  -  A  h  ) 
s 


=  a  g  dt 


1/3  1/2 

Setting  W  =  Ag  h  ,  Equation  119  is  transformed  to 


dW 

W(1  -  W3) 


Referring  to  page  61  of  Grcdshteyn  and  Ryzhik  (1983),  integration  of 
Equation  120  gives 


(120) 


jL 

3 


In 


a28 

~2~ 


t  +  C, 


(121) 


Inserting  initial  conditions  and  transforming  into  the  original  variable 
h  =  H  -  Z  ,  Equation  121  becomes 


H  - 


Z  = 


1/2 


H 

o 


Z 

o 


i2As8 


H 

o 


(122) 


Equations  117  and  122  can  be  used  to  estimate  the  variables  H  =  H(t)  and 
Z  =  Z(t)  . 

Trapezoidal  Breach 


106.  For  trapezoidal  breach,  the  cross-sectional  area  is  defined  as 

A,  =  b(H  -  Z)  +  s(H  -  Z)2  (123) 

b 

where  b  is  the  bottom  width.  During  computations,  the  bottom  width  is 
assumed  as  constant. 

!. i near  case 

107.  Dividing  Equation  62  bv  1  n  64  and  using  Equation  123  yields 


64 


(124) 


i  =  ^rr  [b( h-z)  +  s(h-7)2; 

2  s 


Setting  h  =  H  -  Z  and  separating  the  variables,  Equation  124  becomes 


*  ,  ,,  ,  ,2  a  A 

-a.A  +  bh  +  sh  2  s 
2  s 


Referring  to  page  68  of  Gradshteyn  and  Ryzhik  (1983)  ,  the  integral  of 
Equation  125  is 


(125) 


2  1/2 

b  +  4a^sA 


(b“  +  4a  sA  )^2  -  (b  +  2sh)  . 

In  — r - - — ryx -  =  — —  Z  +  C 

(b^  +  4a  sA  )  '  +  (b  +  2sh)  2s 

2  s 


(126) 


Defining  the  integration  constant  and  inserting  in  Equation  126  results,  after 
some  algebraic  manipulations,  in  the  following: 


\  (  \  1/2  1  T'  2  \l/2 

2s (H  -  Z)  =  <  (  b  +  4a?SA  »  -  b  fb  +  4ct2sAs )  +  b  +  2s(H0  "  ZQ) 


2  \^2 

b  +  4a  SA  -  b  -  2s(H  -  Z  ) 

2  s/  o  o 


2  V^2  (  2  \ 

b  +  4a^ sA  /  +  b  exp  lb  +  4a?sA^  I 


\fb2  +  4a„sA 


Z  -  Z  /  a  A 
o  2  s 


+b+  2s(H  -  Z  )  +  lb  +  4a  sA 


-  b  -  2s(hq  -  Zo)  exp  fb2  +  4a2sAs)  (z  -  W 


(127) 


108.  The  second  equation  required  for  determination  of  variables  1!  and 
Z  can  be  obtained  again  by  subtracting  Equation  64  from  Equation  62. 


-  77  =  -  \r  [k(H  -  7.)]1/2[b(H  -  /.)  +  s  (H  -  Z ) 2 1 

dt  d  t  A 

s 


+  ■*.,[  (}]  -  Z)  j 


6  5 


Writing  Equation  128  in  terms  of  the  hydraulic  head  and  separating  the  vari¬ 
ables  yields 


dh 


h^',‘"(-a  A  +  bh  +  sh2) 
L  s 


s1/2 

-  -  dt 


1 1 2 

By  setting  W  =  h  ' "  ,  Equation  129  transforms  to 


dW 


1/2 


(-a„A  +  bW2  +  sW4) 

/  s 


2A 


dt 


Referring  to  page  67  of  Gradshteyn  and  Ryzhik  (1983) ,  integration  of 
Equation  130  gives 


where 


and 


dW 


A?  +  sW 


dW 


A  +  sW 

O 


1/2 


b"  +  4a„sA 
2  sy 


1/2 


2A 


t  +  C, 


2  2  \ 


1/2 


b“  +  4a„sA 
2  s 


1/2 

A8  =  2  +  l(  b  +  4a2sAs  ) 


109.  Equation  131  can  be  further  simplified  as 


(129) 


(130) 


(131) 


(132) 


(133) 


6o 


where  i  is  the  imaginary  number.  Transforming  back  to  h  and  inserting  the 
value  of  0^  according  to  initial  conditions,  after  some  algebraic  calcula¬ 
tions,  results  in 
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The  system  of  Equations  127  and  135  determines  the  two  unknown  variables  H 
and  Z  . 

Nonlinear  case 

111.  The  nonlinear  case  of  trapezoidal  breach  cannot  be  solved  analyt¬ 
ically  even  for  simple  specialized  coefficients. 


PART  VI:  APPLICATION  AND  RESULTS 


112.  The  performance  of  the  BEED  model  was  tested  for  two  historical  dam- 
failure  cases:  the  man-made  earth-fill  Teton  Dam  in  Idaho,  USA,  and  the 
Huaccoto  natural  dam  on  a  tributary  of  the  Mantaro  River  in  Peru.  The  input 
parameters  for  both  cases  were  taken  from  existing  data. 

Simulation  of  Teton  Dam  Failure 


113.  Teton  Dam,  a  93-m-high  earth-fill  structure,  experienced  failure  on 

5  June  1976.  On  3  June  1976,  leakage  was  detected  at  the  toe  of  the  dam.  By 

early  morning  on  5  June  a  large  leak  caused  by  piping  occurred  AO  m  below  the 

crest  near  the  right  abutment,  and  by  noon  of  that  day  the  crest  of  the  dam 

was  breached.  The  water  in  the  reservoir  was  almost  at  full  capacity  (3.1  x 
3  3 

10  m  ),  and  the  total  mass  of  water  was  released  in  approximately  A  hr,  pro- 

A  3 

ducing  a  maximum  outflow  discharge  of  6.6  x  10  m  /sec.  At  peak  flow,  the 
breach  was  estimated  as  trapezoidal  with  150-m  top  width  and  slopes  1V:0.5H. 
This  failure  event  caused  $70  million  in  property  damage  and  loss  of  six  human 
lives . 

1 1 A .  The  dam  had  a  915-m-long  crest  and  was  composed  of  mixtures  of  clay, 
silt,  sand,  and  rock  fragments  obtained  from  excavations  and  borrow  areas  of 
the  Teton  River  canyon  area.  The  geometric  characteristics  of  the  dam  are 
given  in  Figure  18.  The  reservoir  capacity  and  the  spillway  and  powerhouse- 
outlet  discharge  curves  are  defined  in  Figure  19.  Figure  20  shows  a  map  of 
the  flooded  area,  which  extended  approximately  103  km  downstream  to  near  the 
Shelley  gaging  station  on  the  Snake  River.  The  geometric  and  physical  charac¬ 
teristics  of  various  sections  of  the  flooded  area  are  defined  in  Table  3  (Ray 
and  Kjelstrom  1978). 

Model  calibration 

115.  The  BEED  model  simulation  was  started  after  an  initial  trapezoidal 
breach  with  sides  1V:0.25H  developed  at  the  crest  of  the  dam.  The  slopes  of 
the  dam  were  taken  as  1V:3.08H  and  1V:2.55H  for  the  upstream  and  downstream 
faces,  respectively.  The  initial  hydraulic  head  was  taken  equal  to  1  m  of 
water  flowing  from  the  trapezoidal  breach,  where  the  ratio  of  bottom  width  to 
water  depth  was  0.5.  The  median  soil  particle  diameter  was  taken  as  3  mm,  the 
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Figure  18.  Geometric  characteristics  of  Teton  Dam 


OUTLET  &  POWERHOUSE 
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Figure  19.  Operational  characteristics  of  Teton  Dam 
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Table  3 

Geometric  and  Physical  Characteristics  of 


Cross  Sections  of  the  Flooded  Area  Between 
Teton  Dam  and  Shelley  Gage  on  the  Snake  River 


Distance  from 

Manning ' s 

Station 

Teton  Dam 

Width 

Slope 

Coefficient 

Number* 

km 

m 

% 

of  Friction 

1 

0 

300 

0.0025 

0.032 

2 

7.24 

470 

0.0018 

0.035 

3 

10.46 

4,400 

0.0012 

0.032 

4 

14.00 

4,600 

0.0027 

0.030 

5 

26.39 

10,200 

0.0042 

0.037 

6 

36.36 

5,200 

0.00032 

0.042 

7 

52.45 

5,200 

0.00076 

0.040 

8 

62.59 

3,900 

0.0014 

U.037 

9 

71.76 

1 , 300 

0.0014 

0.035 

10 

75.95 

260 

0.0014 

0.033 

11 

90.75 

140 

0.0014 

0.040 

12 

97.83 

200 

0.0014 

0.035 

13 

103.30 

170 

0.0014 

0.036 

Source : 

Ray 

and  Kjelstrom  1978. 

Station  ' 

Locations  are  shown 

in  Figure  20. 

porosity 

as 

0.2,  and  the  specific 

weight  ol 

f  soil  as  2.5  tons 

(2.3  metric 

tons) /m 

.  Other  soil  characteris 

;ics  were 

assumed  as  follows 

cohesion 

strength 

49, 

000  N/m  ;  angle  of  internal  friction,  40  deg;  and 

Chezy's  coef- 

fieient  < 

1  /° 

of  roughness,  50  m  /sec 

(Ray  and 

Kjelstrom  1978) . 

The  relation  for 

the  Einstein 

-Brown  bed-load  transport  function  for  values  of 

1/’:  higher  than 

2  was  taken  as 

4>  =  139. (136) 


The  coefficient  and  the  exponent  in  Equation  136  were  obtained  by  trial  and 
error,  since  there  are  no  laboratory  or  field  data  in  that  range  of  high  shear 

"7  O 
/  rl 


stresses.  All  data  required  for  flood  roucing  were  taken  from  Table  3.  The 

assumed  inflow  hydrograph  I  =  I  (t)  is  shown  in  the  following  tabulation. 

o  o 


Time  Discharge 

hr  m3/sec 


0 

0.5 

1 

2 


28.32 

84.96 

158.59 

212.40 


3 


192.58 


4 


130.27 


6 

8 

10 

27.5 


67.97 

42.48 

28.32 

28.32 


Simulation  results 

116.  The  computed  outflow  (Figure  21)  indicates  that  the  timing.  shap° , 
and  magnitude  of  the  estimated  outflow  hydroer?rh  compare  quite  well  with  the 
observed  values.  The  higher  peak  outflow,  total  volume,  and  delayed  falling 
limb  of  the  simulated  hydrograph  may  be  attributed  to  the  fact  that  the 
program  continued  the  simulation  until  erosion  reached  the  bottom  of  the  dam, 
while  in  reality,  the  breach  bottom  terminal  elevation  was  about  14  m  above 
ground  level.  At  peak  outflow  discharge,  the  simulated  breach  was 
trapezoidal,  with  a  top  width  of  161  m  and  side  slopes  of  1V:0.675H.  This 
breach  is  very  similar  to  the  field  data  breach  reported  as  having  a  152. 4-m 
top  width  and  side  slopes  of  1V:U.5H.  After  the  peak  flow,  however,  the  model 
estimated  a  total  width  of  295  m,  much  higher  than  a  recorded  width  of 
approximately  200  m. 

117.  Routing  of  the  predicted  outflow  discharge  produced  a  hydrograph  for 
station  13  near  Shelley  that  was  sharper  and  much  higher  than  the  one  recorded 
(Figure  22) .  The  discrepancy  can  be  attributed  to  the  fact  that  the  BEED 
model  used  the  values  of  slope  and  friction  that  were  given  for  the  flood 
channel,  and  not  its  adjacent  floodplain,  which  presumably  had  less  slope  and 
higher  friction.  Another  possible  reason  for  the  inaccuracy  is  the  steepness 
of  the  breach  outflow  discharge,  which  acted  like  a  shock  wave.  Improvement 
can  be  obtained  either  by  using  a  more  accurate  routing  approach,  such  as  the 
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Figure  22.  Outflow  discharge  hydrograph  at  Shelley  station 


St.  Venant  system  of  equations,  or  extending  the  Muskingum  method  to  incorpo¬ 
rate  two-dimensional  flow  behavior. 

118.  To  evaluate  the  dependence  of  the  system  on  the  various  physical 
parameters,  a  sensitivity  analysis  was  conducted  by  changing  the  input 
parameters  within  certain  ranges.  The  basic  characteristics  were  as 
f ol lows : 


r  sn  l/2i 
C,  =  50  m  /sec 

h 

D,_q  =  3  mm 

C  =  49,000  N/m 


4>  =  40  deg 

p  =  0.2 


7  5 


The  relative  difference  in  discharge  (RDD)  is  defined  as 


RDD 


x  100 


(137) 


where  is  the  computed  outflow  discharge  and  is  the  basic  outflow 

discharge  and  both  were  computed  as  presented  in  Table  4. 

119.  Table  4  offers  evidence  that  the  peak  outflow  discharge  is  insensi¬ 
tive  to  most  physical  parameters  except  the  angle  of  internal  friction.  The 
time  of  occurrence  of  peak  discharge  increases  with  increasing  Chezy's  coeffi¬ 
cient  of  friction.  The  opposite  is  true  for  the  median  particle  diameter, 
cohesion,  internal  friction,  and  porosity.  The  side  slope  of  breach  at  the 
time  of  peak  discharge  is  larger  for  low  values  of  cohesion,  internal  fric¬ 
tion,  and/or  porosity.  In  general,  within  a  certain  Hegree  of  accuracy,  the 
model  can  be  considered  as  insensitive  to  variation  of  those  physical 
parameters . 

120.  Special  attention  should  be  given  to  the  determination  of  the  ratio 
between  breach  bottom  width  and  hydraulic  head  and  to  the  expression  for  the 
Einstein-Brown  formula.  Unfortunately,  the  data  are  insufficient  to  define 
these  characteristics  accurately,  and  a  trial-and-crror  procedure  is  required. 

Simulation  of  Failure  of  Huaccoto  Dam 


121.  On  April  25,  1974,  a  landslide  that  occurred  in  Cochacay  Creek,  a 
tributary  of  the  Mantaro  River  in  Peru,  created  the  Huaccoto  natural  dam.  The 
landslide  material  consisted  mostly  of  silty  sand  and  clay  with  D  ^  of  about 
11  mm,  but  there  was  also  material  in  size  up  to  1  m.  The  embankment  was 

170  m  in  height,  and  its  lateral  base  length  was  3,803  m.  The  geometric 

characteristics  of  Huaccoto  Dam  are  represented  in  Figure  23.  The  dam  created 

8  3 

an  artificial  lake  with  a  maximum  capacity  of  8.87  x  10  m  .  The  reservoir 
capacity  curve  is  given  as  Figure  24. 

122.  Huaccoto  Dam  failed  from  overtopping.  During  the  period  6-8  June 
1974,  the  overtopping  flow  created  a  channel  along  the  downstream  face  of  the 
dam  and,  in  the  next  6  to  10  hr,  a  rapid  increase  in  erosion  resulted  in 
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Table  4 

Results  of  Sensitivity  Analysis  for  Teton  Dam 


Parameter 

Value 

Chezy's  coefficient 

of  friction,  m  /sec 

30 

40 

50 

60 

70 

RDD 

5.10 

1.00 

0 

3.60 

3.09 

Time  of  peak  outflow 
discharge,  hr 

0.74 

1.00 

1.26 

1.49 

1.73 

Breach  side  slopes 
at  peak  discharge, 

IV :  sH 

0.675 

0.675 

0.675 

0.675 

0.675 

Median  particle 
diameter,  mm 

1 

3 

5 

RDD 

1.39 

0 

-1.17 

Time  of  peak  outflow 
discharge,  hr 

1.94 

1.26 

1.15 

Breach  side  slopes 
at  peak  discharge 

0 . 67  a 

0 .  o  7  5 

u .  b  /5 

2 

Cohesion  in  N/m 

30,000 

40,000 

49,000 

RDD 

0.48 

5.55 

0 

Time  of  peak  outflow 
discharge,  hr 

1.42 

1.26 

1.26 

Breach  side  slopes 
at  peak  dis. harge 

0.811 

0.675 

0.675 

Soil  internal  friction 
angle ,  deg 

30 

40 

RDD 

18.99 

0 

Time  of  peak  outflow 
discharge,  hr 

1.42 

1.26 

Breach  side  slopes 
at  peak  discharge 

0.811 

0.675 

Porosity 

0.2 

0.3 

0.5 

RDD 

0 

0.33 

-5.18 

Time  of  peak  outflow 
discharge,  hr 

1  .26 

1.04 

0.71 

Breach  side  slopes 
at  peak  discharge 

0.6  75 

0.6  75 

0.5  50 
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a  final  trapezoidal  breach  of  107-m  depth,  200-  to  230-re  top  width,  and  side 
slopes  of  about  1V:1H.  The  peak  outflow  discharge  was  reported  to  be  f roc;  1.0 
:■  104  to  1.80  x  104  m^/sec. 

Model  calibration 

123.  The  BEED  model  started  the  simulation  after  an  initial  trapezoidal 

breach  with  sides  1V:0.25H  developed  at  the  crest  of  the  dam.  The  slopes  of 

the  dam  were  taken  as  1V.-10.56H  aid  1V:12.49H  for  upstream  and  downstream 

faces,  respectively.  The  initial  hydraulic  head  was  taken  equal  to  1  m  of 

water.  The  ratio  of  bottom  width  over  hydraulic  herd  was  0.4.  The  median 

particle  size  was  taken  as  15  mm,  the  porosity  as  0.4,  and  the  specific  weight 

3 

of  soil  as  2.5  tons  (2.3  metric  tons)/m  .  Other  soil  characteristics  were 

2 

assumed  as  follows:  cohesion  strength,  40,900  N/m  ;  angle  of  internal  fric¬ 
tion,  40  deg;  and  Chezv's  coefficient  ->f  friction,  50  m^  ^ / sec  (Ponce  1982, 
Freau  1984).  ihe  Einstein-Rrown  bed-load  function  for  values  of  l/l  higher 
than  2  was  taken  as  constant  and  equal  to  360.  The  inflow  hydrograph  was  con¬ 
sidered  to  be  very  small  and  thus  was  neglected. 

-Simulation  results 

124.  The  computed  outflow  presented  in  Figure  25  shows  that  the  shape, 

timing,  and  magnitude  of  the  computer  l-ydrograpb  resembled  fair!'-  well  the  one 

observed.  Her®  it  should  be  noted  tb°t  higher  simulated  discharge  in  this 

case  does  not  constitute  inaccuracy  for  the  model  because  there  an3  conflict- 

3 

ing  reports  of  the  actual  pea«  discharge  ranging  from  10,000  to  18,000  m  /sec 
with  a  more  likely  value  of  13,500  m' /sec.  At  the  time  of  peak  outflow  dis¬ 
charge,  the  n?.  each  had  a  tcp  width  of  167  m  and  ride  slopes  of  1V:  1.378H. 

Those  results  are  very  close  to  field  reports  that  give  a  top  breach  width  of 
198  m  and  side  slopes  of  The  terminal  breach  obtained  by  the  simula¬ 

tion  was  8 1  m  ».<  :p  and  had  a  bottom  width  of  14.61  m.  The  side  slopes  were 
IV: 1.6711,  indicating  that  mass  failures  of  the  breach  side  slopes  occurred 
following  j . nk  discharge. 

17  4.  'hiring  t  .  1  ihraf  ion  of  the  BEEP  model  t  or  the  lluaccoto  Parr.,  it  was 
noted  tha ’  tin  red- load  function  used  for  I'et-'-n  Pam  (Equation  1  361  was  net 
npplic  Me,  as  t  v  war  giving  very  hi;-!;  erosion  rates  and  ;  ubsequent  1  v  high 
nt  1  fv  discharge.'..  T!  e  effect,  n-  the  various  -.elected  functions  on  the 
magnitude  and  t  in.ii.g-_  ■  outllow  q  i :  ehargi  are  pre-.-nt  i-d  in  t  at  '  e  ’  1 .  wing 
t  .  7:  a  i  i .  1 1  ion. 


Function 


a 


154.55  160.00  171.48  211.12  226.27  320.00 


fm  =  B  1'05  1‘°°  °'9  0,6  0,5  0 


Peak  outflow 

discharge,  m  / sec 

81,000 

78,015 

73,077 

68,021 

54,861 

16,700 

Time  of  peak  outflow 
discharge,  hr 

19.34 

19.56 

19.99 

21.85 

22.80 

27.30 

The  sensitivity  of  the 

model  to  the 

bed-load 

discharge 

formula 

is  an 

indica- 

tion  of  the  need  for  basic  research  on  erosion  processes  under  high  shear 
stresses . 

Analytical  results 

126.  The  uncertainty  associated  with  the  lumped  erosivity  coefficient  a0 
in  the  analytical  solutions  and  the  assumption  of  a  constant  water  surface 
make  the  calibration  of  those  solutions  very  difficult.  Comparing  these  solu¬ 
tions  with  the  governing  equations  of  the  BEED  model,  it  is  evident  that  the 
coefficient  c*„  is  a  function  of  surface  roughness,  cohesion,  angle  of  inter¬ 
nal  friction,  porosity,  specific  gravity,  and  water  viscosity. 

i'll .  To  investigate  the  benavior  of  analytical  solutions,  a  limited  number 
of  tests  were  conducted  for  the  nonlinear  case  of  the  rectangular  breach.  ror 
those  tests,  a  hypothetical  dam  of  the  following  characteristics  was  assumed: 

initial  height  of  dam,  1.  m;  initial  breach  bottom  elevation,  0.90  m;  breach 

2 

width,  0.2  m;  water  surface  area,  1,000  m" ;  and  erosivity  coefficient  equal  to 
0.001  and  0.01  for  the  first  and  second  cases,  respectively.  In  the  first 
case,  after  100  sec  the  breach  bottom  elevation  was  0.975  m  and  the  hydraulic 


head  was  0.026  m.  In  the  second  case ,  after  48  sec,  the  dam  had  eroded  com¬ 
pletely  and  the  hydraulic  head  receded  by  a  maximum  of  0.87  m. 

128.  Changing  the  water  surface  area  from  1,000  to  1O0  m,  setting  the 
breach  width  equal  to  0.1  m,  assuming  an  initial  1  reach  bottom  elevation  of 
0.95  aid  trving  tin-  same  problem  lor  ;  0.1  and  0.01  ,  complete  erosion 


1  !  rVf 

d  in  8  n.in 

i  i ;  the  lit  s t  <  a  e ,  while 

it  took  n  2  min 

in  the 

SOI  I'll 

h » ■ 

peak  hvdr.oil 

i  <•  le  ad  was  0 . 89  in  !  or 

i  ,  -  < 1 .  1  and  ’  * 

.81  lor 

M  1 


129.  Those  few  results  stress  the  fact  that,  before  attempting  to  use  the 


analytical  solutions,  both  the  erosivity  coefficient  and  the  water  sur¬ 

face  area  should  La  properly  calibrated  and  defined.  An  interesting  result 
from  the  analytical  solutions  is  that  peak  outflow  discharge  is  obtained  at 
the  time  the  breach  reached  its  terminal  depth.  This  is  in  agreement  with 
results  from  the  BEED  model. 


PART  VII:  CONCLUSIONS  AND  RECOMMENDATIONS 


130.  The  conclusions  obtained  from  this  study  can  be  summarized  as 
follows: 

a.  With  consideration  of  both  surface  erosion  and  sloughing  of  the 
breach  sides  due  to  instability,  the  BEED  model  satisfactorily 
simulates  the  processes  of  gradual  erosion  of  earth-fill  dams. 

b.  The  accuracy  of  the  BEED  model  is  comparable  to  the  accuracy  of 
other  existing  models,  but  the  physical  basis  of  the  governing 
processes  is  improved. 

£.  The  model  is  not  very  sensitive  to  the  variations  of  the  rough¬ 
ness  coefficient,  the  median  particle  diameter  porosity,  and  soil 
cohesion.  For  changes  within  reasonable  ranges  of  values,  the 
model  produced  a  relative  error  in  peak  outflow  discharge  of  less 
than  5  percent. 

d.  The  internal  angle  of  friction  of  the  soil  is  important  for 
estimation  of  peak  outflow  discharge. 

e.  Low  values  of  porosity,  internal  friction,  and/or  porosity  result 
in  a  wider  terminal  breach. 

f_.  The  ratio  of  bottom  breach  width  over  hydraulic  head  performed 
well  for  values  of  0.5  and  0.4  for  Teton  Dam  and  Huaccoto  Dam, 
respectively.  However,  there  is  no  evidence  that  the  same  values 
can  be  applicable  for  other  cases, 

g.  The  Einstein-Brown  formula  for  values  of  the  dimensionless  shear 
stress  (1/T)  greater  than  2  is  difficult  to  estimate  a  priori, 
and  a  trial-and-error  process  is  required.  The  results  are 
sensitive  to  the  selection  of  this  function. 

h.  The  peak  outflow  discharge  is  obtained  when  the  breach  reaches 
its  terminal  bottom  elevation. 

i.  Sloughing  of  breach  sides  may  occur  even  after  the  occurrence  of 
peak  outflow  discharge. 

j_.  The  routing  part  of  the  BEED  model  overestimates  the  discharge 

hydrograph  .c  downstream  stations.  This  may  be  attributed  to  the 
fact  that  the  breach  outflow  hydrograph  is  very  steep  and  acts 
like  a  shock  wave.  Also,  two-dimensional  effects  of  flow  over 
the  floodplain  might  contribute  to  the  discrepancies. 

131.  Based  on  the  experience  gained  from  this  study,  the  following  recom¬ 
mendations  for  further  study  are  given. 

a.  Extension  of  the  BEED  model  to  include  options  of  other  sediment 
transport  theories,  such  a:;  Schoklitsch,  Mever-Peter  and  Mueller, 
and  the  Bagnold  power  method. 

b.  Application  of  the  BEED  model  to  as  many  historical  cases  as  pos¬ 
sible  for  estimation  of  the  range  of  physical  parameters,  and 
determination  of  the  optimum  technique  for  sediment  transport 

processes . 
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£.  Calibration  of  the  analytical  solutions  according  to  data 
obtained  from  simulation  of  the  BEED  model. 

d.  Extension  of  the  routing  part  of  the  BEED  model  to  incorporate 
two-dimensional  flow  along  the  floodplain. 

£.  Extension  of  the  BEED  model  by  coupling  flood  routing  with 
sediment  routing  along  the  downstream  channel. 
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APPENDIX  A:  USER'S  MANUAL  FOR  BEED-I  MODEL 


1.  The  BEED  model  has  two  versions:  BEED-I  is  written  in  FORTRAN  77 
language  and  is  intended  for  use  in  mainframe  computer  facilities,  while 
BEED-II  is  written  in  ZBASIC  and  is  suitable  for  Zenith  Z100  microcomputers. 
Both  versions  contain  the  same  computational  elements,  and  their  objective  is 
the  same.  However,  version  II  may  impose  memory  storage  limitations,  while 
its  efficiency  measured  in  CPU  time  is  30  to  60  times  less  than  version  I. 

2.  The  program  consists  of  the  MAIN  program  and  six  subroutines: 

a.  MAIN  program  includes  the  iterative  algorithm  for  estimation  of 
water  surface  level  H2  and  breach  bottom  elevation  DZT1.  It  also 
contains  computations  of  the  discharge  QTT  at  the  various  down 
stream  channel  cross  sections  during  flood-routing  processes. 

b.  Subroutine  COVOL  computes  changes  in  reservoir  water  volume 
capacity . 

c.  Subroutine  C0QS01  estimates  spillway  outflow,  powerhouse  outlet 
outflow,  and  inflow  discharge. 

d.  Subroutine  COQBV  estimates  bed-load  sediment  transport  in  volume 
per  unit  width  and  unit  time,  using  the  Einstein-Rrown  formula. 

e.  Subroutine  C0Y2YT  calculates  water  depth  at  the  taiiwater  section 
by  using  the  Newton-Raphson  or  a  fixed-point  iteration  algorithm. 

E_.  Subroutine  C0QB2  estimates  reduced  discharge  coefficients  for 
consideration  of  submerged  flow  conditions. 

£.  Subroutine  COATW  calculates  cross-sectional  area  and  top  width  of 
flood  channel  computational  segments. 


Input  Data 

3.  Input  characteristics  in  both  versions  are  the  same.  Of  course,  the 
FORMAT  statements  of  BEED-I  do  not  apply  to  BEED-II.  For  consistency,  how¬ 
ever,  it  is  suggested  that  the  same  format  structure  be  used  for  both 
versions.  Input  data  are  divided  into  groups.  Each  group  contains  certain 
variables,  as  described  below, 
a .  (’.roup  1  . 

(1)  Physical  characteristics  and  properties  of  dam. 

CM,  DS,  P,  NIU,  PH!,  COH  -  Format  (4F8.A,  F11.7,  2F10.2) 

LH:  Chezy's  coefficient  of  friction  tor  the  breach  in 

172, 

m  /sec 
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DS:  representative  value  of  dam  material  diameter  (usually 

D^q)  in  meters 

SR:  relative  specific  weight  of  dam  material  (y  /y) 

P:  soil  porosity 

2 

NIV:  kinematic  viscosity  of  water  in  m  /sec 

PHI:  angle  of  internal  friction  of  dam  material  in  degrees 

2 

COH:  cohesion  of  dam  material  in  N/m 

(2)  Geometric  characteristics  of  dam. 

LB,  LT,  LC,  HTD ,  HB  -  Format  (5F10.3) 

LB:  longitudinal  length  of  the  base  of  dam  in  meters 

LT:  horizontal  projection  of  upstream  slope  of  dam  in 

meters 

LC:  longitudinal  length  of  dam  crest  in  meters 

HTD:  height  of  dam  crest  in  meters 

HB:  elevation  of  dam  bottom  in  meters 

(3)  Initial  conditions  of  dam-breach  computations. 

HIT,  Z1T ,  Z2MIN,  DTB  -  Format  (4F10.3) 

HIT:  initial  water  level  in  the  reservoir  in  meters 

Z1T:  initial  elevation  of  breach  bottom  in  meters 

Z2MTN:  minimum  permissible  breach  bottom  elevation  (usually 
Z2MIN  =  HB)  in  meters 

DTB:  basic  time  step  for  breach  erosion  computations  in 

seconds 

(4)  Breach  erosion  characteristics. 


L1MIN,  SEC,  EXPO,  X,  ZT  -  Format  (5F10.3) 

LIMIN':  minimum  longitudinal  horizontal  length  of  breach  bottom 

at  the  top  of  dam  in  meters 

SEC:  coefficient  of  sediment  transport  function  (i.e.,  f  = 

I .  /unEXPO 
sec (I /T) 


EXPO:  exponent  of  sediment  transport  function 

X:  initial  ratio  of  breach  bottom  width  to  initial 

hydraulic  head 

ZT :  initial  side  slope  of  breach  section 

15)  Choice  of  iterative  algorithm. 

NO FI  -  Format  (12) 


N0F1 :  If  NOF1  =  0  ,  ^he  Newton-Raphson  algorithm  is  used  for 

computation  of  normal  depths. 

If  N0F1  =  1  ,  the  fixed-point  iteration  is  used  for 
computation  of  normal  depths. 

b .  Group  2 . 

(1)  Computational  time  characteristics. 

TSIM,  DTR  -  Format  (2F11.7) 

TSIM:  time  of  simulation  in  hours 

DTR:  time  interval  for  flood  routing  in  hours 

(2)  Computational  distance  characteristics. 

DXMIN  -  Format  (FI 0.3) 

DXMIN:  maximum  permissible  channel  reach  distance  in  meters 

Note :  If  distance  between  two  consecutive  stations  is 

greater  than  DXMIN,  the  program  generates 
intermediate  sections  until  that  condition  is 
satisfied . 

(3)  Flood  channel  geometric  and  physical  characteristics . 

BV ( I ) ,  ZV(I),  CHV(I),  NMV(I),  SV(1),  D1ST(I),  PRT(I)  - 
Format  (3F10.2,  2F10.5,  F10.2,  F4.1) 

Note :  This  card  is  repeated  until  all  sections  are  input. 

The  last  card  may  say  9999.0  for  BV(I)  and  0.0  for 
the  rest.  The  first  section  is  used  for  the 
tailwater-ef f ects  check  immediately  downstream  of  the 
dam. 

BV(T):  channel  bottom  width  at  section  (I)  in  meters 

2V(I):  channel  side  slope  at  section  (I)  in  m/m  (the  channel 

is  assumed  symmetrical) 

CHV(T):  Chezy's  coefficient  of  friction  at  section  (I)  in 

1/2, 

m  /  sec 

NMV(I):  Manning's  coefficient  of  friction  at  section  (T) 

Note :  Only  one  coefficient  of  friction  can  be  specified. 

The  other  must  be  set  equal  to  O.D.  All  sections 
must  have  the  same  coefficient  of  friction. 

channel  slope  at  section  (I)  in  m/m 

distance  from  dam  to  section  (I)  along  the  channel  in 
meters 

PRT(I):  parameter  to  specify  whether  hydrograph  at  section 

(1)  will  be  printed;  if  yes,  l’KT(l)  =  1;  if  not, 
PRT(I)  =  0.0 
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SV(I)  : 
I)IST(T)  : 


(4)  Flood  channel  initial  conditions. 


QTTLS  -  Format  (F10.2) 

QTTLS:  assumed  initial  discharge  at  the  last  section  in 

m  /sec 

Note :  Initial  discharges  at  the  rest  of  sections  are 

defined  by  interpolation  between  initial  total 
discharge  from  dam  HYD(l)  and  QTTLS. 

c .  Group  3 . 

(1)  Reservoir  characteristics. 

HVL(I),  VOL (I)  -  Format  (F10.2,  F15.2) 

HVL(T):  reservoir  water  level  associated  with  stored  water 

volume  in  meters 

VOL(I):  volume  of  stored  water  when  reservoir  water  level  is 

HVL(I)  in  irt 

Note :  Maximum  number  of  data,  including  last  card,  is  20. 

Last  card  must  say  9999.0  for  HVL(T)  and  010  for 
VOL(I).  First  card  must  correspond  to  dam  bottom 
elevation. 

(2)  Spillway  characteristics. 

HSP(I),  QSP(I)  -  Format  (2F10.2) 

HSP(I):  water  elevation  associated  with  spillway  outflow 
discharge  in  meters 

QSF(I):  spillway  discharge  at  water  elevation  HSP(t)  in 

3, 

m  /sec 

Note :  Same  comments  as  in  (1). 

(3)  Powerhouse  outlet  characteristics. 

H0V(1),  OOV(T)  -  Format  (2F10.2) 

HOV(I):  water  elevation  associated  with  outlet  discharge  in 

meters 

3 

QOV(I):  outlet  discharge  at  water  elevation  HOW.  I)  in  m  /sec 

Note:  3  am  e  comments  as  in  (1). 

C4)  Inflow  hydrograph  characteristics. 

TIKI),  INF(1)  -  Format  (2F10.2) 

TIFfI):  time  associated  with  inflow  discharge  in  meters 

INFCI):  inflow  discharge  at  time  T1F(I)  in  m  /sec 

Maximum  number  of  data,  including  last  card,  is  20. 
last  card  must  sav  9999.0  for  TIFtl)  and  0.0  for 
TNF(T) . 


A  0 


Note: 


Output  Data 


Output  data  are  outlined  below. 

a.  Part  A  -  Basic  Data. 

DTB:  basic  time  step  for  breach  erosion  computation  in 

seconds 

DTR:  time  interval  for  flood  luuLing  in  hours 

SEC:  sediment  transport  equation  coefficient 
EXPO:  sediment  transport  equation  exponent 

l  j  ? 

CH:  Chezy's  friction  coefficient  for  the  dam  in  m  "/sec 

DS :  dam  material  particle  diameter  in  meters 

PHI:  internal  friction  angle  of  dam  material  in  degrees 

COH:  cohesive  strength  of  dam  material  in  X/nT 

ZT:  initial  side  slope  of  breach  section  in  m/m 

X:  initial  ratio  breach  bottom  width  ov°r  initial 

hydraulic  head 

HIT:  initial  reservoir  water  elevation 

Z1T:  initial  breach  bottom  elevation 

DXMIN:  maximum  distance  for  routing  sections 

b.  Part  B  -  Breach  Erosion  Simulation. 

TT :  time  in  hours 

a 

QB2:  breach  outflow  discharge  in  n  /sec 
H  2  ( R  F  L )  :  ratio  (H/-HB)  /  (H.1T-HB) 

Z2  (REL)  :  ratio  (Z2-HB)  /  (.Z1T-HB) 

BRHT:  breach  section  depth  =  HTD-Z2  in  meters 

B2:  breach  bottom  width  in  meters 

B2T :  breach  top  width  -  R2  +  2.0  *  ZT  *  BRUT  in  meters 

c .  Part  C  -  Hydrograph  at  Damsite. 

T:  time  in  hours 

(ITT:  total  discharge  -  QB2  i  OHi’f  •  0012  in  m  '  sec 

d.  Part  P  -  Flood  Pouting. 

T:  time  in  hours 

1.2,1.:  sections  along  thq  ciiannel  reach,  at  which  the 

\ 

hydrogrnph  OTf  (m'/sec''  is  printed 
TT1Y :  lag  time  ot  hydrograph  at  »  neh  sect  ion 


APPENDIX  B:  LISTING  OF  BEED-I  COMPUTER  PROGRAM 


ISN 


ISN 


ISN 
ISN 
ISN 
ISN 
ISN 
ISN 
ISN 
ISN 
ISN 
I  SN 
ISN 
ISN 
I  SN 

ISN 

ISN 

ISN 

ISN 

ISN 

ISN 


ISN 


I 


2 


3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 


1G 

17 

18 

19 

20 
21 

22 


C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


REAL  HVL  (2  0 1  ,VCL<20 )  ,HSP  (  20  )  ,  OS  P  I?  0)  ,  MCU  ! ?0 ) , TCO ( 2 3  >  , 

1  T  IF  (20  t.I  NF(  20  .  IMF  1 ,  INF  2  .'rnvn  (  3  ?  0  3  I  f  , 

2  QINI3000) .QTT( A60.4),BV( 60  )>, ZVI 40  7) ,CHV( 6  V) )  , 

3  NMVI60  0) .  SV< 600)  , Of  ST  (600  ),  (600)  , 

4  QTTPT  I  50)  frjiUfLP,LT,[  C  , K , L  1 , L? , L SSI  ,L  SG2  ,  L  SS7  ,H(>  *  , 

5  LHS  »  NMS  »Kl»NMV2fLHIN*QMA<M  (601  ) ,  TT  P  A  V  (  6 0  )  , 

6  TTTRI600)  ,TTTp(60) 

INTEGER  SUPM.MPP ST ( 5^) 


:*-*<!'*  *  t  * 


program  be 

-  V  EPS  ION  F  OPTRA’ I 


r  n  I 

?r- 


***********  ***  *******************  *******  **!*•************  ** 


************************************** **i 


MAIN  PI  OOP  AM 

**************************************  ***********  ******** 


****************************  *********** 
F  0  P  M  A  T  S 


10  F0RMAT(4F8.4.F11.7,2F1 3.2 > 

15  FORMAT! 5F10.3) 

20  FORMAT! 4F10.3) 

25  FORM  AT  (  5 FI  0 .3  ) 

30  FORMAT! 12) 

35  FORMAT! 2F1 1 . 7  » 

37  FORMAT! F10  .3  ) 

40  FORMAT(3F10.2,2F10.5,F10.?,F4. 1) 

42  FORMAT! F10.2I 
45  FORM AT ( F10  .2  »  F 15  .2  > 

50  FORMAT! 2F1 0.2) 

100  FORMAT!  •  1*  ,4X, ‘DTP  !  S  )  =  •  ,  F  10 . 0 , 5  X ,  *0  TP  (H)  =  ',F11.M 
105  FORMAT! 5 X, »  SEC=« , F7 . 3 , ax , 1  EX PG= • , r 7 . 3 , 5X ,  'CH= ',f7.N 
1  5Xt«DS=’,F8.5) 

110  FORM  AT  (5X,  *PH!=*  ,F8.2,  S^TOItr ',F«.?,5X, '7T=' ,FS.,?, 

I  5X  »  *  X=  * ,F  5 .2  ) 

111  FORMAT!  5X,  'H1T=*  ,F  10.3  ,5X,  *7.1  T*'  ,F  10.3  ) 

112  F  ORM AT ( 5  X . *  OXM  IN= ' »  F  10 .3  ) 

113  FORMAT!  '0*  ,4X, 'RPrACH  FRCSIQN  S I  MOL  AT  I  CM  *  ) 

115  FORMAT! 1X> 

120  FORMAT!  *0*  f3X,  »  TT  ( H ) • , 3 X. ■ 032  !  '*  3  /  S  )  '  ,  2  X ,  '  M  7  ( *'  E  L  )  '  . 

1  3X,  •  Z2  (PEL)  •  ,2X,*  BPIIT  !M)*,4X,'R2  (w|*t3y,'R?T  (  m  >  *  t 
130  FORMAT!  IX,  12F1C.2) 


B1 


I  SN 

23 

145  FORMATI IX, • ZT  NEW=’,F6.7) 

150  FORMATI • P  , IX, ’HYDPOGRAPH  AT  tht  n am  S  It F ’  ) 
155  FORM  AT  I  •  3*  ,  6  (  2  X  ,  ’  T  <H)',2X,'0TT  I  V3/ S)  ?X  )  ) 

ISN 

24 

ISN 

25 

I  SN 

26 

165  FORMATI  6I3X,  F5.2,  r  1 1  .2  ,2X  )  ) 

ISN 

27 

170  FORMATI  • 1*  ,  IX, ’FLOOD  ROUTING'  ) 

I  SN 

28 

175  FORMAT! *0* ,4X, 'T  (H)',15X,'QTT  (^3/S)  at  STA" 

ISN 

2*9 

180  FORMATI  ’O'  , 1  OX  ,  10  (  110,1X1) 

190  FORMAT!  IX, F8.3.4X, 10IF 19. i,ix)» 

***** C  *******************  **  **********  ***** 

DATA  PROCESSING 

***********  ************************■*****>•: 

ISN 

30 

C 

C 

.c 

c 

r 

ISN 

31 

L 

465  RE  AD  I  5,  lOKH.DS,  $f  ,P,N  IU,  PH!  ,r  rn 

I  SN 

32 

READ! 5, 15) LB, L  T,LC  ,HTD,HP 

RE  AD  (5,  20)  HIT,?  IT  ,7  2’M  N,DTP 

ISN 

33 

I  SN 

34 

RE  AD (5,  25IL1MIN.SEC ,EXPC,X, ZT 

ISN 

35 

RE AD(5, 30) NOFL 

ISN 

36 

READ! 5,  35) TSIM.DTP 

ISN 

37 

c 

c 

c 

r 

READ(5,37)DXMIN 

READ  DATA  FOR  riOOD  POUTING  AND  DEFINE 
INTERMEDIATE  S  EOT  IONS 

ISN 

38 

L 

J=1 

ISN 

39 

P  EAD 1 5  •  40 ) BV( J ) *  Zv I J ) »  CHV I J ) , N  MV I J  ) , SV  I  J  ) , D  T  c 
1  PRTIJ) 

ISN 

40 

470  PEAO(5,40IBV2,ZV2,CHV2,NMV2,5V2,ni  ST  ?,  PI’T? 

ISN 

41 

IF(BV2. EQ. 9999.0)00  TO  460 

ISN 

42 

IFIDIST2-DISTI J) .OT.DX’M  N)THrf 

ISN 

43 

NINT=  (DI ST 2- C  1ST! J )  )  /P XMIN+1 

I  SN 

44 

ELSE 

ISN 

45 

NINT=  1 

ISN 

46 

END  IF 

ISN 

47 

J1  =  J 

ISN 

48 

JT=J  +  NI  NT-  1 

ISN 

40 

DO  1000  J=J1,J! 

ISN 

50 

BVIJ  +  l) =BV( Ji  M-IBV2-BVI Jll) *( J+  L-Jl ) /NINT 

ISN 

51 

zv(j  + 1)  =zvi  ji  >  *■  <zv2-z  vi  ji  ) )  *(  j+i- ji  i /r:i  nt 

ISN 

52 

CHVI  J*1  )=f  HV  I  J  1  )  ♦!  CUV  2-f  HV  I  J  1 )  )*(  J  +  l- J  1  )  /’ 

ISN 

53 

NMV I  JU  )=NMV(  JI  )  ♦ (NMV2  -NMV  (  JI  )  )  M  JM-J1  l/f 

ISN 

54 

SV(J+1)  =  5V(J1)t(SV2-SV(  JI)  1  *(  J*  1-  J!  )  /*'!  NT 

ISN 

55 

OIST ( J+l  )  =  DIST( JI  |  +  (0  IST2-DISTI JI  )  ) *( J+l-0 

ISN 

56 

PP.T(  J*l  )=0.Q 

ISN 

57 

1000  CONTINUE 

ISN 

58 

PRT  I  JT+-  l )  =  PRT2 

ISN 

50 

J= J 1+NI NT 

ISN 

60 

GO  TO  470 

ISN 

6  I 

c 

c 

r 

460  NS  - J 

READ  INITIAL  DISCHARGE  AT  LAST  SEE  T 1 0 M 

ISN 

62 

V* 

RCAO!  5,  42)  QTTl.  S 

c 

C  READ  VOLUME  STCPEO  BY  THE  RESERVOIR 
C 


B2 


ISN 

63 

ISN 

64 

ISN 

65 

I  SN 

66 

ISN 

67 

ISN 

68 

ISN 

60 

ISN 

70 

ISN 

71 

ISN 

72 

ISN 

73 

ISN 

74 

ISN 

75 

ISN 

76 

ISN 

77 

I  SN 

78 

ISN 

70 

ISN 

80 

ISN 

81 

ISN 

82 

I  SN 

83 

ISN 

84 

ISN 

85 

ISN 

86 

ISN 

8  / 

c 

c 

c 


c 

c 

c 


c 

c 

c 


c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


1=0 

540  1=1+1 

READ (5, 45)HVL(  I )  , VOL (  I  ) 
IFIHVLI  n.NE.«9oQ.o>Gf;  TP  640 

READ  SPILLWAY  DISCHARGE 


1=0 

560  1=1+1 

RE  AO ( 5  *  50 ) HSP I  I ),QSP< I ) 
IFIHSPI  M  .  NC  .  Qt5<?6. 0  )  GP  TC  560 

READ  OUTLET  DISCHARGE 


1  =  0 

580  1=1+1 

READI5, 50) HOUI I > ,QDU( I ) 

IFIHOUI  I). ME. 9090. 0)00  TP  530 

READ  INFLOW  HVf)ROGPAPH  TO  the  RESERVOIR 

1=0 

600  1=1+1 

READ (5 • 50) T IF l It,  INF(I) 

IF (TIF(I) . NE.^090 .0 )GC  TC  600 

PRINT  HEADINGS 

WRITE(6,100)OTR,DTR 
WRITE  (6 ,10  5  )  S  EC  ,  EXPO.CH.DS 
WPITE(6,110)PHI .COM, ZT.X 
WRITEI6.il  DH1T.Z1T 
WRITE (6 .112  )DXM  IN 
WRITEI6 ,11 3) 

WRITEI6, 115) 

WR(TF(6,»’0) 

WRI1E(6»11 5) 


BREACH  EROSION  SIMULA  TIEN 


*********** *** *******  *******  **********  * 

INITIAL  VALUES 
*************************************** 

TT  =  TIME 

DTR  =  BASIC  TIME  tUTERVAl 

OT  =  VARIABLE  TIHE  INTERVAL 

DZ1  =  BPFACH  EROSION  AT  THE  TOP  PT  THE  f!  A  M  OUPINr  nT 

DZ  2  =  BREACH  E°OSICN  AT  THE  DOWNSTREAM  PACE  DUE  I NC  DT 

OZ T 1  =  BREACH  DEPTH  AT  THE  TOP 

DZT2  =  BREACH  DEPTH  AT  H.F  PCWNSt-FA>'  rsrr 

H2  =  RESERVOIR  LEVEL 


B3 


C 

12 

=  BREACH  BOTTOM  LFVEL 

C 

AA 

=  UPSTPEAM  FACE  ANGLE  WITH  THE  HOT  1 7  0NT  AL 

C 

BB 

=  DOW  NST  RFAM  FACE  ANGLE  WITH  THE  HORIZONTAl 

C 

K 

=  EINSTE IN-BROWN'S  EOUATIC^  COEFFICIENT 

C 

A 

=  DUMMY  VARIABLE  USED  TO  COMPUTE  B2 

C 

B2 

=  BREACH  BOTTOM  WIDTH 

c 

HY0( J I =  OUTFLOW  HYDRGGRAPH 

c 

r 

NT 

=  NUMBER  OF  TIMT  STEPS  USED  FOR  M^nU) 

ISN 

88 

v. 

TT=0. 0 

ISN 

89 

DZ1=0.0 

ISN 

90 

DZ2  =0 .0 

ISN 

91 

H  2=H IT 

ISN 

92 

Z2=  Z IT 

ISN 

93 

DZT1*HTD-Z2 

ISN 

94 

D  Z  T  2=HTD-Z 2 

ISN 

95 

AA=  AT  AN  (  (HTD-HB  I/LT  1 

ISN 

96 

BB=ATAN( (HTO-HR » / ( LB -L T- LC ) ) 

ISN 

97 

K=  36. 0*NIU*NI  U/  (9. 8*DS**  3*1  S*»-  ID 

I  SN 

98 

c 

c 

c 

K=SQRT(2./3.*K)-SQPT(KI 

COMPUTE  SPILLWAY, OUTLET  E  INFLOW  DISCHARGES 

ISN 

99 

CALL  C0QS0KH2.GSP2  ,HSP,  CS  P  ,  QE1U2  ,  HCU  ,  0  CU ,  T  T  ,  INP?  ,T 

1  INF) 

ISN 

100 

QMAX= .0001 

ISN 

101 

A =H2-Z2 

ISN 

102 

c 

B2=X*A 

C  COMPUTE  Q82  AND  CHECK  FOP  TAILWATFP  EFFECTS 
C 

ISN  103  B2S=BVl 1 ) 

ISN  104  ZS^ZVdJ 

ISN  105  CHS=CHVm 

ISN  106  NMS=NMV<1) 

isn  107  ss=svm 

ISN  108  CALL  C0QB2(H2»Z2»B2»ZTtN0Fi»B2S»7S»CHS*NMS,SS,op2, 

1  22MIN, SLBM,Q0U2 ,QSP? ,HBJ 
ISN  109  NT=3600*TSIM/DTfl 

ISN  110  QTT2=QSP2*Qau2*QB2 

ISN  111  HYr(l)=QTT2 

C 

C  ***********  *****************************  ************** 
C  LOOP  TO  CCMPUTC  HYO(J)  -BEGINNING 

C  ****************************************************** 

C 

ISN  112  DO  1482  J=1,NT 

C  *************************************** 

C  VALUES  FPOM  LAST  ITERATION 

C  *♦*  ************************************ 

c 

C  COMPUTE  BRFACH  ROTTC*  WIDTH 
C 

ISN  113  1035  IF (H2-Z2.GT .A )A=H2-Z2 

ISN  115  B2=X*A 


ISN 

116 

IFIQMAX.LT „C82 )QMAX=QB2 

ISN 

118 

D  T=DTB 

ISN 

1  19 

H  1=H2 

ISN 

120 

Z 1  =  Z2 

ISN 

121 

DZT1=DZT1«-DZI 

ISN 

122 

DZT2=DZT2-t-0Z2 

ISN 

123 

CALL  CO VOL (H1»V0L1,HVL,V0L) 

ISN 

124 

QB  1=QB2 

ISN 

125 

QTT1=QTT2 

c 

C  I F  QB2/QMAX.LE. 0.005, DEFINE  QB2  Rv  I NTFR PCL AT  I CN  AM  n  DO 
C  NOT  COMPUTE  DZ 1.0Z2. 

C 

ISN  126  IF(QB1/QMAX.GT.0.005)GC  Tp  1055 

ISN  127  QB2=.0005*QMAX+(QB1-.0005*3MAX)*(360C*TSI**-tT-DT)/ 

1  (3600*TSIM-TT) 

ISN  123  H2=Z2+(QB2/(1.45*B2I )  **  (  2./3.  ) 

r 

C  Ul=WATER  VELOCITY  AT  THE  TOP  PF  THF  DAM 
C 

ISN  129  U1=QB2/  UB2+ZT*  t  H2-Z  2  J)*  ( H2-Z2  I) 

ISN  130  GO  TO  1280 

ISN  131  1055  IFIDZTl .GE .HT0-Z2MIN) THFN 

ISN  132  0ZT1=  HTD-Z2M  IN 

ISN  133  GO  TP  1130 

ISN  134  ENOIF 

C 

C  L 1 =BR E ACM  LENGTH  AT  THE  TrP  OF  THE  PAM 
C 

ISN  135  L  l=LC+0  ZT  1  *  ( 1/T  AN!  AA  )-»l/TAN(BB  )J-0ZT2/$IN(BP) 

ISN  136  IFILl.LT.O.OITHEN 

ISN  137  L  1=L  1MI N 

ISN  138  0ZT1=IL 1-LC+0ZT2/SI N< BP ) )/(  l/TANI  AA »♦ l/TANI  RB) ) 

ISN  139  END  IF 

ISN  140  IFIDZTl .GE .HTD-Z2MIN) THEN 

ISN  141  DZT1  =  HT  D-7.2M  IN 

ISN  142  L 1=0 . 0 

ISN  143  D7.T?  =  S  I  MI  BB  )  *(  LC+DZTi*!  l/TAN(  A  A  ( +  l/TANI  RB)  )  —  L  1  ) 

ISN  144  ENOIF 

C 

C  L2=  BP.  EAC  H  LENGTH  AT  THE  DOWNSTREAM  FACE 
C 

ISN  145  L2=IHTD-HB-DZT1)/$INIBB) 

C 

c  ***********  **************************** 

C  LOOP  TO  COMPUTE  H2  -BEGINNING 

Q  *** *************************  *********** 

C 

ISN  146  1130  H2T=H2 

C 

C  ********************************** ***** 

C  LOOP  TO  COMPUTE  DZl  -BEGINNING 

p  ***$**#♦ *** * ** ************************* 

c 

ISN  147  1140  0Z=DZ1 

ISN  148  Z2=HTD-DZT  1-DZ 1 


B5 


I  SN 

149 

IFIZ2.L  E.Z2MINITHEN 

ISN 

150 

12=12  MIN 

I  SN 

151 

Q  BV1  =  0 .0 

ISN 

152 

0Z1=0 .0 

ISN 

153 

Q«V2=0.0 

ISN 

154 

DZ2=0  .0 

ISN 

155 

GO  TO  3230 

ISN 

156 

C 

c 

r 

ENDIF 

COMPUTE  QB2  ANO  CHECK  FOR  SUBMERGENCE 

ISN 

157 

B2  S=3V ( I ) 

ISN 

158 

ZS=ZV( 1 ) 

ISN 

159 

CHS=CHV  (1) 

ISN 

160 

NMS=NMV(1) 

ISN 

161 

ss=svm 

ISN 

162 

c 

c 

c 

c 

r 

CALL  COQB2(H2,Z2,B2,ZT,NOri,B2S,Z 
1  Z2MIN,  SU8M,Qf)U2  ,QSP2  ,FB) 

**************  *********************** 
DETERMINATION  OF  DZ1 
**********************************  *** 

I  SN 

163 

V- 

C 

C 

r 

Ul=QB2/{ (8  2+ZTMH2-Z2) )*(H2-Z? >) 

COMPUTE  SEDIMENT  DISCHARGE  (VOLJMr  RA 

ISN 

164 

L 

CALL  C0QBV!U1,SR,CH,DS,FSI l,SrC,E 
IF(Z2.NE.Z2MIN»G0  TO  1190 

I  SN 

165 

ISN 

166 

QBV1=0. 0 

ISN 

167 

DZ1=0 .0 

I  SN 

16  8 

GO  TO  1210 

ISN 

169 

1190  DZ1=QBV1*DT/(L1*< 1-P)  ) 

ISN 

170 

IFIDZ1.LE.1 .OIGO  TO  1210 

ISN 

171 

DT=DT/1 0 

ISN 

172 

GO  TO  1190 

ISN 

173 

c 

c 

c 

c 

c 

c 

c 

c 

r 

1210  IF(ABS(CZ1-CZ 1 .GT.0.005IG0  Tn  114 

**********************************  *** 
LOOP  TO  COMPUTE  DZ 1  -END 
***********  ***  *********************** 

****************************  ********* 

COMPUTATION  CF  M2 

************************************* 

ISN 

174 

c 

3230  IFIABSI HI- H2 1 . L T. 0 . 0 1 ) THEN 

ISN 

175 

H5V=  H  1-  2 .3 

ISN 

176 

ELSE 

ISN 

177 

H  S  V=  H  2 

ISN 

178 

c 

c 

r 

END  IF 

COMPUTE  VOLUME  STORED  Bv  THE  RESFPVOII 

ISN 

179 

V-/ 

c 

CALL  COVOL (HSV,VnL2,HVL,V0L ) 

S , OP  2 . 


R6 


C 

r 

COMPUTE  SPILLWAY, OUTLET  C  INFLOW  DISCHARGES 

ISN 

180 

CALL  f  OQSOI  (H2,US«'2,HSP,  CS  P  ,QCU2  ,MOI  ,QCU,  TT+DT  ,  INF  2  , 

r 

1  TIF.INF) 

L. 

c 

r 

COMPUTE  QB2  AND  CHECK  FOR  TAIL  WAT  FP  EFFCCTS 

ISN 

181 

u 

B2 S=8 VI  II 

ISN 

182 

ZS=ZV(  1  ) 

ISN 

183 

chs=chvii) 

iSN 

18  V 

NMS=NMVm 

ISN 

18  5 

SS=SV(  1  ) 

ISN 

186 

CALL  C0CB2 <H2,Z2, B2,ZT ,NOF 1 , 32S , ZS ,C US , NM5 , SS , OP  2 , 

1  Z 2MI N,  SUB M , QQU2 , Q SP2 » HB ) 

ISN 

187 

H2C=H1  +  DT*0.5*(  INFl+INF2-QTTl-QSr'2-QrW2-CR2)*(HSV-HI) 
1  f ( V0L2 -VOL l I 

ISN 

188 

I  F  (  ABS<  H2C-H2J  .LC.  0.005)0,0  TO  1270 

ISN 

189 

H2=H2C 

ISN 

190 

GO  TO  3230 

ISN 

191 

r 

1270  IFIABSIH2-H2T)  .GT.0.01IGC  H33 

c 

**********************************  ***** 

c 

LOOP  TC  COMPUTE  H2  -FND 

c 

r 

************************  *************** 

t 

c 

**************  ************************* 

c 

DETERMINATION  OF  0Z2 

c 

r 

*  **************************************  * 

ISN 

192 

V 

IFISUBM.EQ.IITHEM 

ISN 

19  3 

U2=0 . 0 

ISN 

194 

ELSE 

ISN 

195 

SD=<  HTD-HR I / (L  B-LT -l  C I 

ISN 

196 

r 

Z0=ZT*0ZT1/DZT2 

L 

c 

r 

COMPUTE  Y2  AND  U2  ALONG  THE  DOWNSTREAM  FAC F  or  nr  DAM 

ISN 

197 

V- 

CALL  C0Y2YT  (NOF  l  ,  B2  ,  ZO  ,  CU,  0 . 0  »  SD ,  CB  ?  *  v2  ) 

ISN 

198 

U2=QB2/( <B2*ZD*Y2  >*Y 2  ) 

ISN 

199 

r 

END  IF 

c 

r 

COMPUTE  SEDIMENT  DISCHARGE  (VOLUME  BASIS) 

ISN 

200 

CALL  COCBV(t'2,SP,CM,nS,FSI2,SFr  ,CXPO,0PV2,K) 

I  SN 

201 

IFIZ2.LE.Z2MINITHEN 

ISN 

202 

U2=0. 0 

ISN 

203 

QBV2=0.0 

ISN 

204 

nz2=o.o 

ISN 

205 

ELSE 

ISN 

206 

DZ2=Q  BV2  *DT /  ( L2  *  ( 1  -P  )  ) 

I  SN 

207 

r 

END  IF 

c 

r 

PPINT  RESULTS 

ISN 

208 

1280  TT=TT*Dt 

ISN 

209 

QTT2=0SP2*0OU2*QB2 

B7 


1SN 

210 

TTH=TT/36CC 

ISN 

211 

PH 2= CH2-HB)/! HIT- HR ) 

ISN 

212 

R  Z2= ( Z2-HB ) / ( Z1 T-H3 ) 

ISN 

213 

BRHT=HTD-Z  2 

ISN 

214 

B2T=B2«-2.0*ZT*BPHT 

ISN 

215 

r 

WRITE (6,13  0)TTH  ,QR?,PH?,P72,R0HT,B2,  R2T 

c  ****************************  *********** 

C 

DETERMINATION  OF  LATERAL  SLOPP 

c  ************** **** ******** ************* 

r 

ISN 

216 

L. 

IFIQB2/QMAX.LE.9.O051CC  TO  147} 

ISN 

217 

AP=ATAN( 1/ZTI 

I  SN 

218 

1355 

AP=AP-5/57 .29578 

ISN 

219 

IF(AP)1470,1360,1360 

ISN 

220 

1360 

L  SS1= (HTD- 2 2 )  *(  1/TAN (APJ-7TI 

ISN 

221 

L SS2=(H IT- Z2)*< 1/TAN (AP)-ZT) 

ISN 

222 

LSS3=(H2-Z2)*ZT 

ISN 

223 

H4=  LSS2*(H1T-H2 >*TAN<AP ) / ( H2-Z 2+ (H 1T-M2 ) *  I 1-ZT* 

1  TAN  (  AP  )  ) ) 

ISN 

22  4 

A  1=0.5*11  SS1*IHTD-Z  2)-LSS2*<  H1T-Z2  >> 

ISN 

22  5 

A2=0.5*LSS2*H4 

ISN 

22  6 

A3=0.5*LSS2*(H1T-Z2)-A2 

ISN 

22  7 

A4=0.5*L  SS  3*  (  H  1T-H2-  H4  ) 

ISN 

22  8 

FG=9800*(SR*Al4(SR+PI»42-MSR-(  1-P))*A3-A4) 

ISN 

229 

FH=4900*(I!  1T-H2IMH2-Z2) 

ISN 

230 

NOPM=FG *C0  S ( AP )-FH*S INI  AP ) 

ISN 

231 

IF (NORM.LE . 0.0 INORM=0.0 

ISN 

233 

LHS=FG*SIN  (AP)  ♦FH*COSC  API 

ISN 

234 

RHS= NORM*T AN  IP  HI/57.29  57  8) +C OH *(HTO-Z?)/ SI NIAP1 

ISN 

235 

IF( LHS. LE. RHS )GG  TO  1355 

ISN 

236 

ZT=1/TAN(AP) 

I  SN 

237 

WRITEI6  1 145  IZT 

ISN 

238 

1470 

IF! J*DTB.GT.TT>GC  TO  1035 

ISN 

239 

HYDIJ+1  )-QTTU(  J*DTB-TT4DT)*(QTT?-QTT1)  /CT 

I  SN 

240 

1482 

/- 

CONTINUE 

C  *** 

***********  **************************************  ** 

c 

LOOP  TO  COMPUTE  HYD(J)  -END 

C  ****************************  ******************** 

r- 

ISN 

241 

WRITFI6 ,150) 

ISN 

242 

WR ITE(6 , 15  5 ) 

ISN 

243 

WRITEI6,  U5) 

ISN 

244 

ntt=nt*i 

ISN 

245 

DO  1504  J=1,NTT 

ISN 

246 

THYDI J)={ J-l )*0TB/36O0 

I  SN 

247 

1504 

CONTINUE 

ISN 

248 

WRITE!  6, 165MTHY0I  J),HYD(  J),J=l,NTT) 

C  ********************************************************* 

C 

c  riocD  pouting 

c 

C  ********************************************************* 

c 


B8 


C  PRINT  HEADINGS 
C 

ISN  249  503  WRIT5(6,170) 

ISN  250  WR I TE ( 6 1 17  5  I 

C 

c  **********************  ***************** 

C  INITIAL  CONDITIONS 

C  *************************************** 

C 

C  HYDROGPAPH  IN  SECTION  1 
C 

ISN  251  NP-TSIM/OTP 

ISN  252  1=2 

ISN  253  DO  1550  N= 1  ,NR 

ISN  254  1=1-1 

ISN  255  1545  I F ( I *DTB-N *3 60 0*0 TO .GE. -0.001 ICO  tt  1S4? 

ISN  256  1=1+1 

ISN  257  GC  TO  1545 

ISN  253  1547  QI N IN)  =  HYD ( 1 1  +  ( HYD( I +1 )-HYO( I  )  ) * ( N*3 600* DTP - ( I  li 

1  *OTB)/DTR 

ISN  259  1550  CONTINUE 

C 

C  DEFINITION  OF  THE  SECTIONS  to  BE  '’PINT ED 
C 

ISN  260  M=  1 

I SN  261  DC  1650  J= 1, NS 

ISN  262  IFIPRTI JI.FG.O.OJGO  TC  1650 

ISN  263  MPRST (M)=PRT( J) 

ISN  264  M=  K+l 

ISN  265  1650  CONTINUE 

ISN  266  M T=M- 1 

ISN  267  WRITE  (6  »180  )  (MPRST  (  M  )  ,  M=  1 ,  MT  ) 

ISN  268  WR ITE ( 6  » 11 5 ) 

C 

C  INITIAL  0 IS  CHARGE  AT  EVEPY  SECTION  (LINEAR  I  NT  E°  POL A  T I Cf  1 ) 
C 

ISN  269  DO  1680  J=1,NS 

ISN  270  QTT( Jf2 ) =QTTLS  +  (HYD ( 1  I -CTT LS ) *( NS- J I / (NS- 1 > 

ISN  271  QTT( J+l ,2I=QTTLS+(HYD(l)-QTTLS)*(NS-3-l >/<N5-l) 

ISN  272  QMAXMU  1  =  0.001 

ISN  273  1680  CONTINUE 

C 

c  *********** *** ************************* 

C  LOOP  THROUGH  TSIM  -BEGINNING 

0  ********************************** ***** 

C 

C  AL  PHA1 =KI NE  MAT  I C  WAVE  APPROXIMATION  COEFFICIENT  AT 
C  SECTIONS  (  J  ,  N )  AND  (J,N+1) 

C  ALPHA2  =  KINEMATIC  wave  approx  imat  ion  COEFFICIENT  at 
C  SECTION  (J+1,N) 

C  BETHA=  KINEMATIC  wave  APPROXIMATE  EXPONENT 
C  CA=FLOCD  WAVE  CELEPITV  AT  SECTION  (J,\) 

C  CB  =  FLCOO  WAVE  CELFP I  TV  AT  SECTION  <JfN+lJ 

C  CC  =  FLCOO  WAVE  CELERITY  AT  SECT  II  N  (J  +  lfN) 

C 

l SN  274  00  1950  0=  1, NR 


B9 


I  SN 
ISN 


ISN 
ISN 
ISN 
ISN 
ISN 
ISN 
ISN 
ISN 
ISN 
ISN 
I  SN 


ISN 


ISN 


ISN 

ISN 


ISN 
ISN 
I  SN 
ISN 
I  SN 

ISN 


ISN 
I  SN 
ISN 
I  SN 
ISN 
ISN 
I  SN 
ISN 
ISN 
I  SN 
ISN 
ISN 
l  SN 
ISN 


275 

276 


277 

278 

279 

280 
28  1 
282 

283 

284 

285 

286 
287 

28  8 


28  9 


290 

291 


292 

293 

294 

295 
29  6 


29  7 


299 

300 

301 

302 

303 
306 
306 

30  7 

308 

309 

310 

31  1 
312 
31  3 


QTT ( 1 , 3 ) =0 l N ( N I 
NST=NS- I 

c 

C  LOOP  ALONG  CHANNEL  -BEGINNING 

C 

no  1890  J=  1  i N  S  7 
IFINHVI JI.EO.O.OITHFN 

ALPHA 1=CH VI J  >*SQRT (5 VI J> ) 

ALPHA 2=CHV< J  +  l) *  SORT  IS  VI J+  1)  ) 

BETHA  =  1 .5 
ELSE 

ALPHA  1=SQRT(  SV(  J)  t/NMV<  J  I 
ALPHA2  =  S  OPT ( SV<  J  +  l )  * /NWV ( J*l) 

BETHA-S. / 3 . 

END  IF 

CA=B£THA*ALPHAl**  (  1  .0/  RETHA)  3  (  ;tt(  ),2I/!W(JM 
1  BETHA-1  .OI/BFTHA) 

CB=BETHA*Al.  PHA  I**  I  1 . 0/ BE  THA  )  *  (  or  T  (  J,  3  )  /R  V  (  J  )  ) 
l  **(< BETHA-1 .OI/BETHA) 

CC.  =  BETHA*ALPHA2**f  1.0/RETHA)*<0tt(  J*  1  ,  ?  I  /  *»V  «  J  +  1  )  I 
1  **<  I  BCTHA-1 .OI/BETHA  ) 

C 

C  COMPUTE  FtOCO  WAVE  CEIEPITY  AMO  UNIT  WIDTH  0  T  S r  H  A  r  G  c 
C 

C= (CA+C  B+C  C  )  /  3  .0 

QN  = ( QTT ( J, 2 ) /BVf  J ) +OTT( J , 3 ) /RVI J) + QTT ( j+ 1 , ?  > / 

1  BV  I  J-*-  1 )  I  /  3 . 0 
C 

C  COMPUTE  PARAMETERS  FOP  mijskINCUM  ME  THrH 
C 

IF (C .EQ.O. OITHEN 

QTT(J-*-l»3)-=QTT(  J+1,2) 

ELSE 

K1  =  (D I ST( J+l ) -DI ST( J  J) /C 

X=0.5*(1.0-2*QM/( ( SV (  J ) +  S V ( J  *1  I  )*C  *( n[ sT( J+  I ) 

1  -OISTfJim 

IF ( X. L  T. 0. 0) X  =  9. 0 

C 

C  COMPUTE  COEFFICIENTS 
C 

C4=36C0*0TP  /K  1+2 . 3*1  1- X  ) 

Cl  =  <3  600*OTR/K1+2.0*X)/C4 
C  2=(3600*DTR/Kl-2.  C*X>  /C4 
C  3=  12  .0*  (1  -X  )-360  3«OTP  /K  ll/C  4 

QTTI  J*1  ,  3)  =C  1*QTT(  J,  2>+C2*QTT<  J  ,  3)  *C3*  y  t  ( JM  ,?  ) 
IFIQMAXM  <  J*1 ) .LT.QTTI J* 1 ,3)t  0M\XM( J+ 1) =^TT(J*1  ,3  I 
IF(QTT(J*1,3I .LT.0.0 (THEN 
P PINT*  , • J=»  ,J 
PP I  NT*, 'PI S  T ( J  1= • .PI  ST  t  J) 

PRI  NT*  »,QTT(J«-1,3)  =  *  , QTT I J ♦ 1  ,3  ) 

Pp I  NT* , *  MODIFY  SPACF  AMD/OP  TImr  INtfpv'Lc' 

STD  P 
E NO  IF 
ENpir 


BIO 


ISN 


314 


1890  CONTINUE 


I  SN 
ISN 
ISN 
1  SN 
ISN 
ISN 
I  SN 
ISN 
ISN 
ISN 
ISN 
ISN 
ISN 
ISN 
I  SN 
ISN 
I  SN 

ISN 
ISN 
I  SN 
ISN 
ISN 
ISN 
ISN 
ISN 
ISN 
ISN 
I  SN 
ISN 
ISN 
I  SN 
ISN 
ISN 
I  SN 


ISN 


315 

316 

317 

318 

319 

320 

321 

322 

323 

324 

325 

326 

32  7 

328 

329 

330 

331 

332 

333 

334 

33  5 
336 

33  7 

338 

339 

340 

341 

342 

34  3 

344 

345 

346 
34  7 
34  8 


349 


C 

C 

c 

c 

c 

c 

c 

c 


LOOP  ALONG  CHANNEL 


C  TR  Av=  PEA  N  WAVE  CEL  ER  I  TV 
TTP  AV=ACC  UMULA  TED  TINT  OF 


•crin 


r  *  *  *  *  * 


TtMvfL  "r  thp  wAVr 


TR=N*DT  R 
M  =  1 

DO  1920  J= 1 » MS 

IF(  PRTI  J) .EC. 0.0  )G0  TO  192^ 

QTTPT (Mi=GTT( J, 2) 

M=  M+-1 

19?0  CONTINUE 
MT=M-l 

W  RITE  (6, 190  ITP,  (QTTPT  (M)  ,  3=1,*''-; 

DC  1940  J=1,NS 

QTTI J, 1 )=QTTI j,2) 

QTT I J  ,2  >  =QTT (  J  ,  3  ) 

1940  CONTINUE 
1990  CONTINUE 

T  TTP,  ( 1  )  =0 . 0 
DO  2100  J= 1 » NS T 

CTP  A V=R  E  THA*  (  SQP  T(  S  V(  J))/N‘,V(JI)**(l.C/PFTHAK 
1  (0.5*CMAXM(J )/BV(J»  )**( (RETHA-1  .0  )/prTHA> 

TTPAV(J*ll-(DrST(J+lJ-OIST(J'  l/CTPAV 
TTTR  (  J*  1 1 -0  .0 
2100  CONTINUE 

DO  2300  J=l,NST 
DO  2400  1=1 . J 

TTTP  ( J*1  i  =TTTR(J«-1  )  ♦TT^AVI  1*1) 

2400  CONTINUE 
2300  CONTINUE 

PRINT*,  ’ - LAC  T  I  pr S - * 

M=1 

00  2500  J=1,NS 

IE (PRT( j| .EC. 0.0 (GO  TC  2900 
TTTP(M)=TTTP( JI 
M=  W+l 

2500  CONT I  Aft/E 
MT=M-1 

WRITEI6 ,190)TR, (tTTP(m)  i ,mt ) 


C 

C 

C 

c 

C 

c 

c 

c 

c 

c 

c 


STOP 


END  OF  t  H  r  A  I  N 

,t*t**«**tt*tn* 


prop  \  •* 


f  t  *  t  t 


I 


SUB  ROUT  INE  COVCL < HSV , V L,  HV  L  »  V  ri  L  ) 


C  ********* ***********  ****** ************* 
C  DETERMINATION  OF  RESERVOIR  VOL'INE 
0  ************************** ************* 
c 


ISN  2  REAL  HVL ( 2 0 1 t VOL ( 2^ ) 

ISN  3  1=1 

ISN  4  IF  (HSV.LT.  HVL(  nJTHFN 

ISN  5  VL=0 . 0 

ISN  6  ELSE 

ISN  7  3050  IFCHSV.LE.MVH  1*1  I  )GC  TO  306T 

ISN  8  1=1+1 

ISN  9  GO  TO  3050 

ISN  10  3060  VL=VOL(  I»  +  (VOL(  I  +  l  )-Vf  L(  m*CHSV-MVL  (  r  I) /(MVL  (I  +  1) 

l  -HVL(I)) 

ISN  11  ENCiF 

ISN  12  PETURN 

ISN  13  END 


812 


1 


SUBROUTINE  CQQSOI I H2 ,CSPl ,HS°, CSP, COLT  ,  unj,  OOU  ,TS  T 

IINFL  |T  I  F  » I  NF  ) 


ISN 


I  SN 


C 

C 

C 

C 

c 

c 


************************************  **  * 

ESTIMATION  OF  SP I LL WA v t DUTL E T  £ 

INFLOW  DISCHARGES 


REAL  HSPI2  0)  ,QSP(2D)  >HC‘U(  21")  ,QnU(20) 

I  I  NFL 


ISN 

3 

I  SN 

4 

ISN 

5 

ISN 

6 

ISN 

7 

3090 

ISN 

8 

ISN 

9 

I  SN 

10 

3100 

ISN 

11 

C 

C  COMI 

c 

ISN 

12 

ISN 

13 

ISN 

14 

ISN 

15 

ISN 

16 

3130 

ISN 

17 

ISN 

18 

ISN 

19 

3140 

ISN 

20 

C 

C  COMI 
C 

ISN 

21 

ISN 

22 

ISN 

23 

ISN 

24 

ISN 

25 

31  70 

ISN 

26 

ISN 

27 

1  SN 

28 

3180 

ISN 

29 

ISN 

30 

COMPUTE  SPILLWAY  DISCHARGE 


1=1 

IFIH2  .L  T  .MSP  I  I  )  ITHEN 
QSPL=  0.0 
ELSE 

IFIH2 .LE.HSPI I  +1 > )0D  TO  3100 
1=1  +  1 

GO  TO  3090 

QSPL=QSPU)  +  <QSP( I  +  l I-QSPI I ) )*(  H2- 
i  -Hspnn 
ENDIF 


1=1 

IF  IH2.L  T.HOUI  1 1  )  THEN 
QO UT=  0  . 0 
ELSE 

IFIH2  .LE.HOUI  1  +  1  l)f.O  TC  3140 
1=1+1 

GO  TO  3130 

QOUT=QOU  (!)♦<  QOUt  I  +  l  I-COUU  )  )  MH2- 
L  -HOU (I) ) 

ENOIF 


1=1 

IFITSI.LT.TIFI I  I ITHEN 
1  NF  L=  1  NF  (1 ) 

ELSE 

IFITSI.LE.TIFI  I  +  IUGO  TO  3130 
1  =  1+1 

GO  TO  3170 

I  NFL*  INF  1  n+IINFU  +  1 1-IMFI  II  I*  (T  SI 
l  /ITIFI  I  +  D-TIFIIH 
ENDIF 
RETURN 


! 


,  T  [F  {  20  )  ,  I  *'F  (  20  1  , 


HSP  1 1  I  I/MSPII  ♦  U 


nr:u( » I  I  /(hpij<  I  +  1  I 


-T  IF  I  I  I  I 


ISN  1  SUBROUTINE  CC'QRVI  II,  SF  ,  f  U ,  OS  ,  F  S  I  ,  SEE  ,  E  X  PO  ,QF  V  ,  K  ) 

C  4c  4444******  4c4c4t*4c#*j|4c*4c*4c*4i4c4:^*<>-k*4;4*44i» 

C  COMPUTATION  OF  SEDIMENT  DISCHARGE 

C  (VOLUME  BASIS)  fiv  USING 

C  E INSTEIN-BROWN' S  EQUATION 

C  **  +**  *  *****  **#  ****#*4e*4c*4r4t4r*  ******  *4,*** 


C 


ISN 

2 

REAL  K 

ISN 

3 

fs  i=u*u/(<  sr  -1  )*r.n*c  h*ds  ) 

I  SN 

4 

IFIFS  I  .LE.0.047 )  THEN 

ISN 

5 

FHI =0 . 0 

ISN 

6 

ELSE  IFIFS  I  .L  E  .  0.0562  )THEfl 
FHI=14*FSI-0.1R8)**1.4 

I  SN 

7 

ISN 

8 

ELSE  IF (FS  l.LE  .2 . 3) TMFN 

ISN 

9 

FH1=40.0*FS  1**3 

ISN 

10 

ELSE 

ISN 

11 

FHI  =  SEf *F S I **E  XPO 

ISN 

12 

ENDIE 

ISN 

13 

Q8V=FHI *K* ((SP-1)*9.8*DS**3)**0.5 

ISN 

14 

RETURN 

I  SN 

15 

END 

b!4 


1 


SUBPOUT  INE  COY2YT( N0F1  f  B  2$  ,  7  S  ,  CHS  r  NM  S  ,  3  3  ,  Q  S  ,  v?c  I 


r 


I  SN 


I  SN 

2 

ISN 

3 

ISN 

4 

I  SN 

5 

ISN 

6 

I  SN 

7 

ISN 

8 

ISN 

9 

ISN 

11 

I  SN 

12 

ISN 

13 

I  SN 

14 

ISN 

15 

ISN 

16 

I  SN 

17 

ISN 

18 

I  SN 

20 

ISN 

21 

I  SN 

22 

ISN 

23 

ISN 

24 

ISN 

25 

r 

C  ************************ ********** ****  * 

C  COMPUTATION  OF  NORMAL  DEPTH 

C  ***********  *** *************************** 

c 

PEAL  NMS 

IF (NOF 1 . EQ . 1IOC  TO  3571 
C 

C  NEWTQN-RAPHSGN  NETHCO 
C 

IF ( NMS .EC. O.OITHEN 

Y2S=(  »QS/fCHS*B?Sn**2/SS>**(  1  ./3.  I 
ELSE 

Y2S=(  1.44*NMS*NMS*QS*QS/  I B ’S** 2. *7*SS )  )** ( 1 .  /  3.  ) 
ENDIF 

3532  IF  (NMS.NE.O.O)CHS=(  < R2 S+ ZS ♦V3S ) * v? V  (R2S*2*V2S* 

1  (  l *ZS  *  *2 I  **D .5 ) )**(  1./6.I/NMS 

FY=QS*SQPT < B2S*2*Y2S*SCPT ( 1 +7S  **2  I  )-OUS*<fB?S*7S*v25) 
I  *Y2S)  **1.  5*SQP  T(  SS) 

r  PY=  QS*  SQP.T  (  (  I  +  2S**2)/(P?o  +  ?*jQPt(  l*  Z  S  **  2 )  *  v2  5  )  )  - 
1  i. 5*C  H  S* SQP  T  C S  S  *  I B2S  +  Z  S  * v  2  S • *  Y  2  S )  * ( 82S+?*7S*V?S ) 

Y  2C= Y2S-FY/FP Y 

IF  (ABSI  Y2C-Y2S  )  .LF  .0.035  >GO  rr  aaeo 
Y2S=Y2f 
GO  TO  3  532 
C 

C  FIXED  POINT  ITERATION  MrT MOD 
C 

35  71  Y2  S=  0 .0  01 

3572  IF(NMS.NE.O.O)CHS=  (  (B2S+  ZS’tv2>  I  *  y?  S/  (  P  ?S  ♦?  *  v2  S  * 

1  (1*7.S**2  )  **0.5  I  >**(  1  ,/6.  I /NMS 

Y2C=OQS*iJS*(Q2S+2*y2S*(  1  +  Z  5  *  *2  )  *  *0 . 5  )/<CHS*niS* 

1  <B2S*7S*Y2SI**3*SS)  )**<  1./3.  ) 

[F(ABS(  Y2C-Y2S  I.l  E.  O.OT5  )GC  TO  5583 
Y2S=Y2C 
GO  TO  3  5  72 
3580  RETURN 
END 


1 


ISN 


ISN 
ISN 
I  SN 
ISN 
ISN 
I  SN 
ISN 
ISN 


ISN 
ISN 
ISN 
ISN 
ISN 
ISN 
ISN 
I  SN 


2 

3 

4 
6 

7 

8 
9 

1C 


1 1 
12 

13 

14 

15 

16 

17 

18 


SUBPOUT  INC  CC0B2IH2 ,72,82  ,  7T,NOF  1 , 82  S,  7 S  ,C.HS  ,N>'S  ,  SS, 
1  QB2, 72MIN,SMBM,aOU2,0SP2,H8 > 

r.  *********  ,*****«*»********‘**i«*»it**»t 
C  COMPUTATION  OF  QB2  AND  TAILWA^FR 
C  EFFECTS  DUO  TO  SUBMERGE NCC 

C  *********************************-****** 

c 

REAL  NMS 
INTEGER  SUBM 
IFIZ2.LE.Z2MIN) Z2=Z2MIN 
IF (H2-Z2.GT.O.OJGO  TC  3410 
QB2=  0.0 
GO  TO  3418 

3410  QB2=1.45*B2*(H2-Z2)**1 .5*1 . 1 5* 7 T*  (M2 -Z?) **2. 5 
QS=QB2*GOU2*QSP2 
C 

C  CHECK  FOP  SURMERGENCT 
C 

CALL  C0Y2YT(N0r  l,82S,7S,rHS, N‘1S,  SS,Q  S,  YT) 
IF{YT*H8-Z2.LE.0.67*JH2-Z2I)Gf  341?) 

QB2=QR2*( 1 -27. 8* ( !  yt*HR-Z? ) / I M2  -  72 ) -0 . 6^ )** > ) 

SUBM= 1 
GO  TO  3420 
3418  SUBM-0 
3420  RETURN 
ENO 
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APPENDIX  C:  LISTING  OF  BED-II  MICROCOMPUTER  PUG  CRAM 


10  ' 

20  '  t  ************************* *1  *******  **************  t.****  ********  *  t**  ttttt*  *  * 

30  ' 

40  '  P  R  G  G  R  A  11  BRED  IT 

5  0  •' 

60  ' 

70  '  •  v  E  r.-  6  i  n  n  r  p  a  s  I  c  - 

80  ' 


90  ' 
1 00 
1  10 
120 
1 30 
l  40 
150 
160 
1T 
130 
1  90 
200 


************************************************************************** 

M  A  I  M  P  O'  G  8  9  o  ri 


*********************************** 


************************************** 
BATA  PROCESSING 
************************************** 


210  OPTION  BASE  1 

220>  DIN  MOL..  (20  •.  ,  VOL  '20'  ,  MSP  <20'  ,  QSF-  <  20'  ,  HCU  '  2*'"  .  QfMJ  (  2'‘"  ,  tif  '20'  ,  IMF  '20’  ,  Tuvp  <  j  on 
O'  ,  HYD  '  3  000'  .ClTN’tOOO'  ,  QTT  '  300 . 2  '  .  B’V  '  2'X' '  ,  ;’V  '  2  O'-’ '  .  CH1.1 ' 200 '  .NMV'201'"  ,  SO  <  COO  1  ,DIST 
(200'  ,  RPT  '  200 '  .  F’RST  (20'  .  QTTC  T  .'20'  .  DMA  <  m  ■  2'V. ;  .  ttrav  <  2rt0 '  ,TTTP'20'"  ttjp  i  200  •• 

220  READ  CH,  PS,  SR,  P,  NIU,  FT  .  CGH 
240  READ  LB.LT, LC.HTP.HB 
250  READ  HIT. 21T, Z2NTN.DTB 
260  READ  L ININ. SEC. EVP0.X.2T 
270  READ  NOR  1 
280  READ  TSIN.DTP 
290  READ  P*MIN 
300  ’ 


310  '  READ  DATA  EOF'  FLOOD  ROLLING  AND  PE  R  1 NE  I  NT'RNEP  I  ATE  SECMONS 


740 

r.'EAD 

EV 

f  1 

■  * 7 

11  <  ' 

,  CH’.-' '  J  '•  .  NN7  ( J>  .  c v  /  T  '  ,  DIRT  /  .*  '  ,  r  R  T  '  J 

—  r~  - 

PEAL' 

EV 

r « 

2\*2 

.  CHV 

2.NM07.SV2, PISr2. RRT2 

760 

IF  B 

OD¬ 

QO 

C'O 

THEN 

500 

T  p  D 

IST 

2 

1*  7  2 

7  (  J  1 

DVM1N  T’-CM  NIMT-ir.’T  i  1  r  ICT2-D  I  E'r  z 

7g 

7  1  --.J 

Qi'i 

j  t  - ; 

*NI 

NT 

-  1 

4 

rnc 

’-..I 

« 

T  n 

JT 

41  ‘ 

r-o  4 

■<  1  ' 

-  tl 

V  i  1 

i  \  +  r 

R'.'2-E '  1 1  ’  ■  *  ,1  +  1  ’  1  NIT 

4  7>- 

j 

+  } 

-  - 

’7  '  .1 

I  >  +  ' 

2V2-ZV  1  -  '  *  '  J  *■  1  -  J  t  '  'MINT 

4  7'* 

CH’: ' 

.7  i  *. 

CHV 

'll' 

•*  C.H','2  -CHV'  .11  ’  '  *  'J>  1-  J  1  n:rr 

4  4  ' 

NNV  ' 

- 

NMV 

'  ’  1 

+  •  NMVC-NM7  Ml  '•  '  *  .  J+t  1  1  ■  ’MINT 

470 

■~7  !  .1 

4  i  • 

— 

7  '  ,7 

;  ’  '  • 

S'.'2-EV  '  2  ’  ’  *  '  1  -2  1  -  'MINT 

460 

D:  37 

•  1  + 

1  , 

-  L*  I 

-  T  •  2 

1  >  +  -  r  IRT2  DIS'r  Ml  '  '  *  M  +  l  -21  '  TINT 

4"i 

F-PT  . 

7  +  i 

<\  • 

49i> 

INF  v  t 

J 

4 

C'.c  3  / 

,7  7  * 

:  > 

-F  F 

T2 :  J 

-21  TINT:  GOTO  TK>’> 

50'‘>  MR  = J 


Cl 


OJ  cr 


5 1 c 

520  '  PEA V  T  N  l  r  T  Ml  0  i  EC  MAFGE  AT  OAT  T  TEC  T  I  Of  J 

54/';  PEAT  OTTL=- 
550  ' 

5. AO  ’  PECC  VOLUME  SLOPED  FA  tHl  r  EC f:P’.,r:  1  - 


50:')  PpAD  m'.'L  ■  I N  VCL  f  T  ;  T  r  H‘ 
600  ' 

61'*  ’  PEAT  CT'  ILL  D I  EC. H  AF  Up 


640  KEAD  HEP  M  >  ,  02* '  I  '  :  I F  *J-:T  ■  1 
6^0  1 

A.-'*'  ’  p PAT  CUTLET  DICCHAPGC 


"Hfri  a  1 


7  10 


^EAD  Hni 1  •  1 \  CCU  >  1  •  :  IF  hpm  8  3  ooor  "''v-ir*-'  ’  • 
'  PE  AD  IMP ’.OH  H,  TFCCFhF  w  -v;:  THE  r-  E'CEF  vn  I P 


,  c.pT  r‘:  •  .1 .' 


-CTf" 


1  =  1 

KT-r  TJF  '  2  .  ,  j'jr 
■  pcTtjT  HE  AC  TUP.’ 


1;  T  T  ' 


"HEM  :  '  ; : 7CT7 


SO  7 

nt 

877 


EC/ 

e.-v  ■ 


-  r.’1 

ro  —  ' ' 

'’4  7 

C/Cr 


3T' 


1 1'  1  ■ 


!.F'7  IMT:  I.  PRINT 

t_RP'  71"  "PIP  f  5  ]  =  "  ;  DTP.  :  Lc'c  I N"  "7_rr  r  • 1 1  •■  "  ;  V  T-  : !  f  r-  T*  IT 

:.cp.  r  rj  "  ’  "EC  -  "  ;  BE'  .  :  LC'I:  r  N  7  "7  >  r  -  "  ;  F  M'-p,  :  ;  .FFPC  "PH-  "  :  CM,  :  "..PPIN"  "7S=  "  :  PE  :  1_C'PI 


"0H--  "  5  crus,  :  l  PRINT 
1  T  :  Lc'v!mt 


i  fp::r  "phi  jpmi  , .  i.  pr  irj-" 

•_  Pr  I  NT  "HI"- "  ;  HI  T,  :  "  ’  IT^' 

LF'F  I  NT  "  CTf  T  M-  "  ;  CTHItJ:  rr 
LC"RIMT  "PREACH  EPQCJDN  7  INULA"  1  7M"  :  '_.r'FI,'(T 
I  RR  T*  IT 

L PRINT  "  ”  "U’  PP7  '"IT  7’  H7  :  C '  7  7 

!_F  r  j  t  |T 


;  7  '  .  :  nc  y  >  |T  "  >  =:  '•  ;  v  ;  >_  r  T  I  f  1’ 


p.p-HT  I  HI  FO  : 


-  V  ■  -  -  *>  r  ^ 


BREACH  EROSION  SIMULATION 

***************  Mtiiuuimti 

************************************ 

:  M  ;  '  !  A  A  '  "  F  7 

**********  *************  ************* 


•  P  7-7  17  !'  t  ME  I  r 1  "E R'.'Al 

•  VAPtAr:  r  t  t me  rrrRVA. 


C2 


i  -  r  ’T; 

1050  '  re-: 

lOoj  '  H“ 

i.'  ’  • 

\  r~-’  i  ' 

[  O-'T  '  r  p 
111'1  '  A 

n  "  '  '  h  •  r 

;  j  nr  '  r.'-r 

1  iC'  ' 

l  1  .“0  HG-Hl T 
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1  1  H>  l  >  *tl 


-Breach  eeqehjm  <r  rMr;  the  or  -hp  i:  ;  •  i  r- 
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-:fj'.lMBFF  QF  '' [MR  'RP-  ■  1 — H I.  L-'  -  . 


•' H rr*-HB  '  I  r  A  FT  AIR  u[[,  A 
■ !  ;  pet rc:  •  *  .■  sr  ;  : -  ■  . 


1  ot-’R  i ;Tt.  '-f  T  i  :  wa  .  ,  O'  T  1  ft  '.inflow  LT  SC  har-SEE 

17?  T  =.  t  t  : SCOUR  4‘F'  i'_  ;  .0'  r=OFR  i.  :  Q.JLV  --CfjL,"r :  1  r  JF  E  =  T  rJF- !. 
1  ? A  '  -  -Hr  ;•[  :  =  -  *A 


1  •  GGtTF'JIE  QR7  AML-  r’U.r.D  ‘*CR  ’iV1  Wh’’>  LF  F  u  - "  - 

-  ?  —  f[s 


1  C-CVOB  AT"-'  •• 

j  7<: fjr  =  v/t  ’.-.ORATE  7M  T-TB  '  :  CH1  7- QRR  7  *0’?'  'J  +GF-.  :  m  • !  -  0 i  ’ 

-  7->0  * 

)  -  [  0  ’  M************* ********************************** 

ito  -  !.?nt"  r?  tomf-citf:  h-t>  :  •  -beginning 

)•*?.;  •  tutttttl  w*  WHttt*  Ht  I  Hill  UH  tit  HUM  Mt 

1  74<->  • 
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I'?.'-  '  '  'A!..'  Ip  R  FROM  l  AST  T  TFF-ArIilN 
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********************* ************** 

1  400 

]  4  )  0 

-OMR'.irr  RF f  AC  H  RCTTCN  WII'-th 

1  470 

J  4~0 

’  IF 

H7-  77  A  THEM  A=M?-rO 

1  440 
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QM AT  OB 7  -HEM  QMA-«-0B7 

1  4 

DT 

-O'B 

1  4”’- 

:  H[ 

■  HP:  7  [  ■-  7  7:  r  7T  [  --  r  7  T  ;  +  77  •  7  7-7-7  7 1 1  -*  777;  >!'■. 

TOOVvI  ■!■ 

t  ■  :  .  L'Ll  ”  '.’1 

T  ~r  t 

o r  r 

? 

1  4?r 

1  4"'~ 

:r  op:  omat  rproji.  mb.  t--.  imu-frr 
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■'MEM  ON  AT  ."MG' THEN  1‘V  .<  ■ 

C3 


1520  QB2=. 0005*QMAX+ < QB 1 - . 0005*QMAX ) * ( 3600 *TS IM-TT-DT) / < 3600*TS IM-TT ) 

1530  H2=Z2+ (0B2/ ( 1 . 45*B2) ) '  (2/3) 

1540  ’ 

1550  '  U1=WATER  VELOCITY  AT  THE  TOP'  OF  THE  DAM 
1560  ’ 

1570  U1=0B2/ ( (B2+ZT* (H2-Z2) ) * (H2-Z2) ) : GOTO  2470 
1580  IP  DZT1  HTD— Z2MIN  THEN  DZT 1 =HTD-Z2M IN: GOTO  1740 
1590  ’ 

1600  ’  L1=BREACH  LENGTH  AT  THE  TOP  OF  THE  DAM 
1610  ’ 

1620  L1=LC+DZT1 * < 1 /TAN (AA) +1 /TAN (BB) ) -DZT2/SIN (BB) 

1630  IF  LI  0  THEN  L 1 =L 1 M I N : DZT 1 = < L 1 -LC+DZT2/S I N ( BB ) ) / ( 1 / TAN ( AA ) + 1 / TAN ( BB ) ) 

1640  IF  DZT1  HTD-Z2MIN  THEN  DZT 1 =HTD-Z2M IN: L 1 =0: DZT2  =  S I N ( BB ) * ( LC  +  DZT 1  * ( 1 / TAN < AA ) 
+1 /TAN (BB) ) -LI ) 

1650  ’ 

1660  '  L2=BREACH  LENGTH  AT  THE  DOWNSTREAM  FACE 
1670  ’ 

1680  L2= (HTD-HB-DZT 1 ) /SIN (BB) 

1690  ’ 

1700  ’  ******************************************** 

1710  '  LOOF  TO  COMPUTE  H2  -BEGINNING 

1720  '  ******************************************** 

1730  ' 

1740  H2T=H2 
1750  ' 

1760  ’  ******************************************** 

1770  ’  LOOP  TO  COMPUTE  DZ1  -BEGINNING 

1780  ■’  ******************************************** 

1790  ’ 

IS 00  DZ=DZ1:Z2=HTD-DZT1-DZ1: IF  Z2  =Z2MIN  THEN  Z2=Z2MIN: GOTO  2070 
1810  ’ 

1820  -  COMPUTE  QB2  AND  CHECH  FOR  SUBMERGENCE 
1830  ’ 

1840  GOSUB  4890 
1850  ’ 

1 860  '  ************************************** 

1870  •  DETERMINATION  OF  DZ1 

1880  '  ************************************** 

1390  ’ 

1900  U 1 =QB2/ ( (B2+ZT* (H2-Z2) ) * (H2-Z2) ) 

19  10  ’ 

1920  •’  COMPUTE  SEDIMENT  DISCHARGE  /VOLUME  BASIS) 

1930  ’ 

1940  U=U1:  GOSUB  4  7‘,0 

1950  FSI 1 =FSI : IF  Z2=Z2MIN  THEN  QBV1 =0: DZ 1 =0: GOTO  1990 
I960  QBV 1 =QBV 

1970  DZ1=QBV1*DT/ (LI* ( 1-P) ) 

1980  IF  DZ1  l1  THEN  DT=DT  '10: GOTO  1970 
1990  IF  ABS  <  DZ 1 -DZ )  .005  THEN  1800 

20|-|0  ’ 

2010  '  ************************************** 

2020  '  COMPUTATION  OF  H2 

2030  ’  ************************************** 


2040  ’ 

2050  '  COMPUTE  VOLUME  STORED  BY  THE  RESERVOIR 
2060  ’ 

2070  IF  ABSIH1-H2)  .01  THEN  HSV=Hl-2'  ELSE  HSV=H2 
2080  GOSUB  4400: VQL2=VL 
2090  ' 

2100  '  COMPUTE  SPILLWAY. OUTLET  ?<  INFLOW  DISCHARGES 
2110  ’ 

2120  TSI=TT  +  DT:  GOSUB  4500:  QSP2=QSF'L :  Q0U2=Q0UT :  I  NF2=  I  NFL 
2130  ’ 

2140  ’  COMPUTE  QB2  AND  CHECH  FOR  TAILWATER  EFFECTS 
2150  ' 

2160  GOSUB  4890 

2170  H2C=H1+DT* . 5* ( INF  1  +  INF2-QTT 1 -QSP2-Q0U2-QB2) * (HSV-Hl > / (V0L2-V0L1 ) 
2180  IF  ABS (H2C-H2) >. 005  THEN  H2=H2C:G0T0  2070 
2190  IF  ABS (H2-H2T) >.01  THEN  1740 
2200  ’ 

2210  ’  ************************************** 

2220  '  LOOP  TO  COMPUTE  H2  -END 

2230  '  ************************************** 

2240  ’ 

2250  '  ************************************** 

2260  '  DETERMINATION  OF  DZ2 

2270  ’  ************************************** 

2280  ’ 

2290  IF  SUBM= 1 1  THEN  U2=0:G0T0  2370 

2300  B2S=B2: CHS=CH:NMS=0 1 : SS=  ( HTD-HB ) / (LB-LT-LC) 

2310  ZS=ZT1*DZT1  'DZT2:0S=QB2 
2320  ’ 

2330  '  COMPUTE  Y2  AND  U2  ALONG  THE  DOWSTREAM  FACE  OF  THE  DAM 
2340  ’ 

2350  GOSUB  5060: Y2=Y2C 
2360  U2=0B2/ ( <B2+ZS*Y2) *Y2) 

2370  U=U2 
2380  ’ 

2390  ‘  COMPUTE  SEDIMENT  DISCHARGE  (VOLUME  BASIS) 

2400  ' 

2410  GOSUB  4770 

2420  FSI2=FSI : IF  Z2=Z2MIN  THEN  U2=0: QBV2=0: DZ2=0: GOTO  2470 
2430  QBV2=0BV : DZ2=QBV2*DT / (L2* ( 1  — P > ) 

2440  ’ 

2450  ’  PRINT  RESULTS 
2460  ' 

2  ".70  i  I  =TT  +  DT 

2480  QTT2=QSP2+aOU2+OB2 

2490  TTH=TT/3600 

2500  RH2= (H2-HB) ' (H1T-HB) 

2510  RZ2= (Z2-HB) / (Z1T-HB) 

2520  BRHT =HTD-Z2 
2530  B2T =B2+2 ' * ZT *BRHT 

2540  LPPINT  USING  "  ” ; TTH; : LF'RINT  USING  "  ######. #  " j  0B2: RH2: RZ 

B2T 
2550  ' 


2 : BRHT ; B 


C5 


2560  ’  ************************************** 

2570  '  DETERMINATION  OF  LATERAL  SLOPE 
2580  '  ************************************** 

2590  ' 

2600  IF  QB2/QMAX  =.005  THEN  2740 
2610  Ap=ATN  (l/ZT) 

2620  AP=AP-5/57. 29578 
2670  IF  AP  O  THEN  2740 

2640  LSS1= (HTD-Z2) *  < 1 /TAN ( AP) -ZT ) : LSS2= (HIT-22) * ( 1 /TAN ( AP) -2T) : LSS7= ‘ H2-Z2! *ZT 
2650  H4=L332* (HI T-H2) *TAN ( AP) / (H2-Z2+ (HI T-H2' * ( 1 -Z  T * T AN ( AP )  )  ) 

2660  Al  =  . 5* (LSS1*  < HTD-Z2 ) -LSS2* (H1T-Z2) ) : A2= . 5*LSS2*H4 : A3=. 5*LSS2* (H1T-Z2' -A2 
2670  A4=. 5*L3S7* (H1T-H2-H4) 

2680  FG=9800* (3F:*A1  +  f SR+P ) *A2  + (SR- ( 1 -P!  ) *A7-A4!  : FH=4900* (H1T-H2) * (H2-Z2> 

2690  NORM=FG*COS <AF) -FH*SIN (AP) : IF  NORM  =0  THEN  NORM-O 

2700  LHS=FG*SIN ( AP) +FH*COS < AP)  : RHS=NORM*T AN (PHI /57 . 29578 ) +COH* ( H T D - Z 2 )  /SIN ( AP) 
2710  IF  LHS  =RHS  THEN  2620 
2720  ZT= 1 /TAN ( AP ) 

2730  LPRINT  "ZT  NEW=";ZT 
2740  IF  J*DTB  TT  THEN  1430 

2750  HYD ( J+l ) =QT  T 1  +  ( J  *D  IB-T  T  +  DT ) * ( 0TT2-QTT 1 > /DT 
2760  NEXT  2 
2770  ’ 

2780  ’  **** *********************************** 

2790  ’  LOOP  TO  COMPUTE  HYD(J>  -END 
2800  '  *************************************** 

2810  ’ 

2820  LPP I  NT: LPRINT 

2830  LPRINT  "HYDROGRAPH  AT  THE  DAM  SITE” 

2840  LPRINT  "  T  [HI  OTT  [M3/5I  T  [H3  QTT  [M3/SI  T  [H]  OTT  [M3. S3  T  [HI 
QTT  [M3  'SI" 

2850  LPRINT 

2860  FOR  J  =  1  TO  NT  + 1 

2870  THYD  ( J )  =  ( J- 1 ) *DTB/7600 

2880  LPRINT  USING  "######.##  THYD ( J ); HYD ( J ) ; 

2890  NEV  T  2 
2900  ’ 

2910  '  ************************************************************************* 

2920  - 

2930  •  flood  routing 

2940  ' 

2950  ’  ************************************************************************* 

2960  ’ 

2970  '  PRINT  HEADINGS 
2980  • 

2990  LPRINT: LPRINT: LPRINT 

3000  LPRINT  "  FLOOD  ROUT  I NG " : LPR I  NT 

3010  LPRINT  "  r  [HI  OTT  [M7/S3  AT  STATION  #  ".‘LPRINT 

3020  ’ 

3030  ’  ***************************************** 

7040  ’  INITIAL  CONDITIONS 

3050  ’  ***************************************** 
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3060  ’ 

3070  '  HYDROGRAPH  IN  SECTION  1 
3080  ' 

3090  NR=INT (TSIM/DTR) 

3100  1=2 

3110  FOR  N= 1  TO  NR 
3120  1=1-1 

3130  IF  I*DTB-N*3600*DTR  =-. 001  THEN  3150 
3140  1=1+1: GOTO  3130 

.150  Q I  N  ( N )  =HYD  (I)+(HYD(I+1)-HYD(I>)*  <  N*  *  DT  R-  (  I  -  1  )  tDTB'  /DTE 

3160  NEXT  N 
3170  ' 

3180  ’  DEFINITION  OF  THE  SECTIONS  TO  BE  PRINTED 
3190  ’ 

3200  M= 1 : FOR  J=1  TO  NS 
3210  IF  PRT  ( J  >  =0  THEN  3230 
3220  FRST (M) =PRT ( J ) : M=M+1 
3230  NEXT  J 

3240  LPRINT  "  - 

3250  MT=M-1 

3260  FOR  M= 1  TO  MT 

TC^O  LF'R  I  NT  USING  "  ##  PRST(M); 

3280  NEXT  fi 
3290  LPRINT 
3300  ’ 

3310  '  INITIAL  DISCHARGE  AT  EVERY  SECTION  (BY  LINEAR  INTERPOLATION) 
3320  ' 

3330  FOR  J=1  TO  NS: QTT < J, 2) =QTTLS+ (HYD ' 1 > -QTTLS) * (NS-J) / (NS-1 ) 

3340  QTT (J+l , 2)=QTTLS+ (HYD ( 1 1 -QTTLS) * (NS-J-1 ) / (NS-1 > 

3350  QMAXM < J )=. 001 : NEXT  J 
3360  ’ 

3370  ’  ***************  * ********* * * *** **********  * ** *  *  * 

3380  '  LOOP  THROUGH  TSIM  -BEGINNING 

3390  '  ********************************************** 

3400  ’ 

3410  ’  ALPHA 1=P I  NEMATIC  WAVE  APPROXIMATION  COEFFICIENT  AT 
3420  -  SECTIONS  (J.N)  AND  <J,N+1> 

3430  '  ALPHA2=I INEMATIC  WAVE  APPROXIMATION  COEFFICIENT  AT 
3440  -  SECTION  (J+l.N) 

3.450  ■  BETHA  =1  I  NEMATIC  WAVE  APPROXIMATION  EXPONENT 

3460  ’  CA  =FLOOD  WAVE  CELERITY  AT  SECTION  (J.N) 

3470  '  CB  =FLOOD  WAVE  CELERITY  AT  SECTION  rJ.N+1) 

3480  ’  CC  =Fl.DOD  WAVE  CELERITY  AT  SECTION  (J+l.N) 

3400  ’ 

3500  FOR  N= 1  TO  NR 
3510  OTT ( 1 , 3) =QIN (N) 

3520  ’ 

3530  '  ********************************************** 

3540  '  LOOP  ALONG  CHANNEL  -BEGINNING 

355(3  ’  ********************************************** 

3560 

3570  FOP  J=1  TO  NS-1 
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3580  IF  NMV ( J ) =0  THEN  ALFHA1 =CHV ( J ) *SQRT ( SV ( J > ) : ALPHA2=CHV ( J+ 1 ) *SQRT ( SV ( J  + 1 > ) : BE 
THA= 1.5: GOTO  3600 

3590  ALPHA  1 =SQRT (SV ( J)  > /NMV ( J)  : ALPHA2=SQRT (SV< J  +  l )  ) /NMV  < J  +  l > : BETHA=5/3 
3600  CA=BETHA* ALPHA  1  ( 1 /BETHA) * (QTT ( J , 2) /BV ( J ) )  '  ( (BETHA-1 ) /BETHA) 

3610  CB=BETHA*ALFHA1 '  (  1 /BETHA)  *  ( QTT  <  J , 3  > / BV ( J ) )  (  (BETHA-1  )  /  BETHA) 

3620  CC=BETHA* ALPHA2 " ( 1 /BETHA) * (QTT ( J+l , 2) /BV ( J+ 1 > )  ( (BETHA-1 ) /BETHA) 

3630  ’ 

3640  "  COMPUTE  FLOOD  WAVE  CELERITY  AND  UNIT  WIDTH  DISCHARGE 
3650  ’ 

3660  C=(CA+CB+CC) /3 1 

3670  QM=  (Q  ( J ,  2)  /BV  ( J  )  +Q  ( J  ,  3)  /BV  ( J  )  +Q  ( J  +  l ,  2)  /BV  ( J  +  l )  )  ,’Z  ' 

3680  ’ 

3690  ’  COMPUTE  PARAMETERS  FOR  MUSKINGUM  METHOD 
3700  ’ 

3710  IF  C=0  THEN  QTT ( J+l , 3) =QTT < J+l , 2) : GOTO  3850 
3720  Kl= (DIST (J  +  l ) -DIST  (J) ) /C 

3730  X  =  . 5*  < 1-2 ' *QM/ (  ( SV ( J ) +SV ( J+ 1 ) ) *C* (DIST (J  +  l ) -DIST (J ))) ) 

3740  IF  X  O'  THEN  X=0‘ 

3750  ’ 

3760  ’  COMPUTE  COEFFICIENTS 
3770  ’ 

3780  C4=3600*DTR/Kl+2 1  * ( 1-X) 

3790  Cl= <3600*DTR/Kl+2 1  *X  > /C4 
3800  C2= ( 3600*DTR/K 1-2 ' *X) /C4 
3810  C3= (2 ' * ( 1-X) -3600*DTR/K1 ) /C4 

3820  QTT  r J  +  l , 3) =C1*QTT ( J , 2) +C2*QTT ( J, 3) +C3*QTT ( J  +  l . 2' 

3830  IF  QMAXM ( J+l ) <QTT (J+l , 3)  THEN  QMAXM ( J+l ) =QTT ( J+l , 3) 

-840  IF  QTT (J  +  l, 3) <0*  THEN  LPRINT  " J  =  " ; J : LPR I  NT  "DIST (J )="; DIST <J >: LPRINT  "0TT i J 
<■1,3)=";  QTT  ( J+l  ,  3)  tLPRINT  "MODIDY  DXMIN  AND/OR  DTR":STOP 
3850  NEXT  J 
3860  ’ 

3870  ‘  ***************************************** 

3880  ’  LOOP  ALONG  CHANNEL  -END 

3890  ’  ***************************************** 

3900  ’ 

3910  ’  CTRAV=MEAN  WAVE  CELERITY 

3920  •  TTR A V= ACCUMULATED  TIME  OF  TRAVEL  OF  THE  WAVE 
3930  ’ 

3940  TR=N*DTR 
3950  M= 1 

3960  FOR  J=1  TO  NS 
3970  IF  PRT(J)=0  THEN  3990 
3980  QTTPT (M)=QTT(J, 3) :M=M+1 
3990  NEXT  J 

4000  LPRINT  USING  "  ###.##  " ; TR; 

4010  FOR  M- 1  TO  MT : LPR I  NT  USING  "  #«####.«  QTTPT (M) NEXT  M: LPR I  NT 

4020  FOR  J=1  TO  NS: QTT ( J, 1 ) =QTT ( J, 2) : QTT ( J, =QTT ( J, 3) : NEXT  J 

4030  NEXT  N 

4040  TTTR ( 1 ) =0 1 

4050  FOR  J=1  TO  NS-1 

4060  CTRAV= BETHA* ( SORT ( SV ( J )  ) /NMV  ^  J ) )  < l  1  'BETHA) * ( . 5* QMAXM ( J ) /BV <  J )  '  ((BETHA-1  B 
ETHA) 
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4070 

4080 

4090 

4  1  00 

4110 

4120 

4  1  30 

4  1  40 

4  1  50 

4160 

4170 

4  1  80 

4190 

4200 

4210 

4220 

4230 

4240 

4250 

4260 

4270 

4280 

4290 

4300 

4310 

4320 

4  330 

4340 

4350 

4360 

4370 

4380 

4’90 

4400 

4410 

4420 

4430 

4440 

4450 

446* 1 

447  0 

4480 

4400 

4500 

4510 

4520 

4530 

4540 

4550 

4560 

4570 

4580 


TTRAV (J  +  l ) = (DIST ( J+l ) -DIST ( J) ' /CTRAV 
TTTR (J  +  l ) =0 1 : NEXT  J 
FOR  J=1  TO  NS- 1 

FOR  1  =  1  TO  J: TTTR (J+l ' =TTTR (J  +  l> +TTRAV( 1  +  1)  :  NEXT  I:NEXT  J 

LPRINT  " - LAG  TIMES  - " 

M=1 

FOR  J=1  TO  NS 
IF  RRT < J ) =0 1  THEN  4160 
TTTR (Mi =TTTR ( J ) : M=M+1 
NEXT  J 

LPRINT  USING  ”  ###.##  ";TR; 

FOR  1*1=  1  TO  MT:  LPRINT  USING  "  ######.# 


;  TTTF (M)  ;  : NEXT  M;  ; LPRINT 


**:M*#******#**#***#**»****************** 

LOOP  THROUGH  T3IM  -END 
***************************************** 

STOP 

************************************************************  ************* 
END  OF  THE  MAIN  PROGRAM 

************************************************************************* 

END 

************************************************************************ 

S  U  B  R  0  U  T  I  N  E  5 


************************************** 

DETERMINATION  OF  RESERVOIR  VOLUME 

************************************** 


1  =  1 

IF  HSV  HVL  '  I  >  THEN  VL  =  0:  GOTO  448'1 
IF  HSV  HVL'I+1'  THEN  1=1+1: GOTO  4460 

VI.  =  VOL  '!'  +  'VOL  <  1+  1  ■  -v'UL  'l'i«  'HSV-HVL  '  I  ' 
RETURN 


DETERMINATION  OF  SP I LL WA v . OU TLE T 
AND  INFLOW  DISCHARGES 


1  COMPUTE  SR1  ILL  WA  i'  DISCHARGE 


1  =  1 

IF  H2  HSPrn  THEN  QSP  L  =0 :  GOTO  4fc4.': 


;HVL  I  +  1  '  -  HVL  ■'  1 
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4590  IF  H2>HSF'(I+1>  THEN  1=1+1. -GOTO  4590 

4600  QSPL=QSP  (  I  )  !  QSF'  (1  +  1)  -QSP  (  I  )  )  *  ( H2-HSF  (  I  )  )  /  (HSF  <  I  +  1  )  -HSF  '  I  )  ) 

46  1  0  ' 

4620  '  COMPUTE  OUTLET  DISCHARGE 
4630  ' 

4640  1=1 

4650  IF  H2  HOU(I)  HEN  QOUT=0:GOTO  4^  10 
4660  IF  H2:H0U'I+1)  THEN  1=1+1: GOTO  4660 

4670  QOUT  =  QOL ( I )  +  (QOU ( 1  +  1 > -QOU '1)1* ( H2-HQU (Ill  /  (HOU ( I  +  1 ) -HOU  < I )  ) 

4680  ’ 

4690  •  COMPUTE  INFLOW  DISCHARGE 
4700  " 

4710  1=1 

4720  IF  TSI'TIF(I)  THEN  I NFL= I NF ( II : GOTO  4750 
4770  IF  TSI  TIF'I+n  THEN  1=1+1: GOTO  4730 

4740  I NFL= I NF  ( I )  +  <  I NF  < I  +  1 ) -  I NF ( I )  >  * ( TS I -T 1 F  (  I)  )  .  '  T  I F  '  I  +  1 ' - r I F ( I)  ) 

4750  RETURN 
4760  ' 

4770  ’  ************************************** 

4780  ’  COMPUTATION  OF  SEDIMENT  DISCHARGE 

4790  ’  (VOLUME  BASIS)  BY  USING 

4800  '  EINSTEIN-BROWN’S  EQUATION 

4810  ’  ***** ********************************* 

4820  ’ 

4870  FSI=U  2/ ( (SR-1 > *CH  2*DS) : IF  FS1  =.047  THEN  FHI=0':G0T0  4866 

4840  IF  F3I  =.0562  THEN  FH T = ( 4*FSI - .  1 88 >  '  1 . 5: GOTO  4060 

4850  IF  FSI  =2'  THEN  FHI=40*FSI  7  ELSE  FHI=SEC*FSI  EXPO 

4860  QBV=F  H  I  *)  *  <  <SF;-1  )  *9. 9*DS  7>  '  .  5 

4870  RETURN 

4880  ’ 

4890  ’  ************************************** 

4900  ’  COMPUTATION  OF  QB2  AND  TAILWATER 

4910  ’  EFFECTS  DUE  TO  SUBMERGENCE 

4920  ’  ************************************** 

4970  ’ 

4940  IF  72  =  2 2M I N  THEN  Z2  =  Z2MIN 

4950  IF  H2-Z2  -O'  THEN  QB2=0 '  : GOT  0  5070 

4960  QB2=1 . 45*B2* < H2-Z2)  1 . 5+ 1 . 1 5* ZT* < H2- Z 2 )  2.5 

4  7  70  ’ 

4980  ’  CHECt  FOP  SUBMERGENCE 
4990  ' 

5000  B2S=BV l|): 7S=ZV 1 1 ) : NMS=NMV < 1 ) : SS=SV ( 1 ) : 0S=0B2+Q0U2+CSP2: CHS=CHV II': NMS=NMV < 
1  ) 

5010  GOSIJB  5060:  YT=Y2C 

5020  IE  YT  +  HB-Z2  .b"t  (H2-Z2)  THEN  QB2=QB2*  1  1  -2~ .  8*  (  (YT  +  HB-Z2)  ■'  '"H2-Z2)  67'  7  '  :  ?U 

BM= 1 ' : GOTO  5040 
5070  9IJBM=0  1 
5040  RETURN 

r  1 6 1 ;  '  ************************************** 

5070  •  COMPUTATION  OP  Y2.YT 


5090  If  NOF 1 = 1  THEN  5070 
51 'TO  ' 

5110  ’  NEWTON— RAF'HSON  METHOD 
5  1  00  ' 

5170  IE  NMS=0  THEN  YDS  =  (  (  QS  /  (  CHS* BOS  1  )  O/SS)  (1/7-1  ELSE  YDS = ( 1 . 44  * NMS  '  0 * OS  0  '  ( BO 
S  0.67*  SS l  >  (1/7) 

5140  IF  NMS  O'  THEN  CH5= ( ( BDS+ ZS* YOS ) * YOS/ ( BDS+D* YDS* ( 1 + ZS  0 '  .  5 )  )  '  (  1  / 6 1 / NMS 

515'/  FY=QS* (BD3+D*YDS* < 1  +  Z3  0>  ,5>  . 5 -CHS* (  ( BDS  +  ZS* YDS '  *  YOS >  1.5*35  .5 

5160  FPY=QS* ( 1 + ZS  0)  . 5/ ( BOS+O* ( 1 + ZS  0)  .5*Y0S)  . 5- 1 . 5*CHS* ( SS* ( BDS+ZS* YOS ' * YOS ) 

. 5* <BDS+D*ZS*YDS) *SS  .5 
517 0  YOC=YOS  FV-FFY 

5180  IF  ABS ( YOC -YDS )  .005  THEN  YOS=YOC:GOTO  5140 

5190  GOTO  5070 
5000  ' 

5010  '  FIXED  POINT  ITERATION  METHOD 
5070  YOS =.001 

5040  IF  NMS  O'  THEN  CHS= (( BDS  +  ZS* YDS )* YDS/ < BDS+D* YDS* < 1 + ZS  ' 0 5 )>'< 1 /fe WNMS 
5050  YDC=(QS  0* ( BOS+O*  YOS*  < 1 +ZS  D>  . 5) / (CHS '  0* ( BDS  + ZS* YDS '  3*SS1)  (1/7) 

5060  IP  ABS ( YOC -YOS )  .005  THEN  YDS=YDC:GOTO  5040 

50 7 o  RETURN 
5080  ' 


57-00 

5710 

5700 

5770 


************************************** 

DATA 

************************************** 


5740 

DATA 

50  *  0. 

0070 ,2.5,9 

.  00,  IE-6,  40,  49000 

5750 

DATA 

577. 4 

090 . 0 ,  1 0 .  <! 

. 1605. 19, 1530.03 

5760 

DATA 

1617.95 

1610.95. 

1570. 07, 300 

5770 

DATA 

1.0,1 79 

99,  i .  o,  <: 

.5,0. 05 

578'/ 

DATA 

1 

5790 

DATA 

T-6,  0. 1 

>377777 

5400 

DATA 

1 000 ■ * 

» 

54 1  0 

DATA 

30* 

'.  < 

>.  l . 

0.0,  0.072, 

0.  < 

>005, 

o , 

1 

54  70 

DATA 

4  7( 

)  ,  / 

,  o , 

0.0.  0.075, 

0.  ( 

>018, 

7240, 

5470 

DATA 

440* 

>.  ' 

i  ,  0  , 

0.0,  0.070, 

0.  t 

>010. 

10458. 

7, 

54  40 

DATA 

4y‘>i 

>  .  < 

s  o , 

0.0,  0.070, 

0 .  ( 

>027 . 

14807, 

4 

5450 

DATA 

1  020i 

>.  ' 

i,  ?■;, 

0.0,  0.077, 

o .  < 

’042, 

06388. 

o 

5460 

DATA 

520f 

5  ' 

>,  0, 

0.0.  0.040, 

7. OE-4, 

76767 , 

6 

5470 

DATA 

520< 

) .  » 

1 ,  0 , 

0.0,  0.040, 

7. 6E-4, 

50457, 

o 

5480 

DATA 

790' 

) ,  r 

> .  o . 

0.0,  0.077, 

0 .  c 

>0  1  4  , 

60590, 

8 

5490 

DATA 

1  3  Of 

> .  f 

» ,  0 , 

0.0,  0.075. 

0 .  ( 

>01  4, 

71761 . 

o 

5500 

DATA 

06' 

) .  * 

• ,  0 , 

0.0,  0.077, 

o  < 

>0  1  4 , 

75945, 

10 

55  1  0 

DATA 

14. 

i ,  l 

* ,  0 , 

0.0,  0.040, 

0 . 00 1 4 , 

90748, 

o 

5500 

DAT  A 

30* 

)  .  ( 

1  ,  o , 

o.O,  0.075, 

0 . oo 1 4 , 

97807, 

0 

55 

DATA 

1  7 

)  „ 

) ,  1*1 , 

0.0.  0.075, 

0 . 00 1 4 , 

107098, 

1  7 

5540 

DATA 

0999 

,  0 , 

0 . 0 ,  0 . 0 

o .  0  , 

o , 

o 

5550 

DATA 

00 .  o 

5560 

DATA 

1570. 07 

0 

557m 

DATA 

1 579 , 04 

0. 0E6 

Cl  1 


5580  DATA  1546.86,  11.01E6 

5590  DATA  1560.58,  38.33E6 
5600  DATA  1569.72,  61.68E6 
5610  DATA  1577.34,  88.11E6 
5620  DATA  1584.96,  123.35E6 

5630  DATA  1591.06,  160.80E6 

5640  DATA  1598.68,  202.65E6 
5650  DATA  1606.30,  251.UE6 
5660  DATA  1622.76,  397.27E6 
5670  DATA  9999,  O 
5680  DATA  1616.96,  0 
5690  DATA  1618.71,  991.2 
5700  DATA  1621.54,  5097.60 
5710  DATA  1625.19,  5097.60 
5720  DATA  9999,  0 
5730  DATA  1566.45,  0 
5740  DATA  1569.72,  566.40 
5750  DATA  1579.52,  718.11 
5760  DATA  1601.29,  849.60 
5770  DATA  1620.88,  960.86 
5780  DATA  1622.76,  960.86 
5790  DATA  9999,  0 


5800 

DATA 

0 , 

28.32 

5810 

DATA 

1800, 

34.96 

5820 

DATA 

3600 , 

158. 59 

5830 

DATA 

7200. 

212. 40 

5840 

DATA 

1 0800 , 

192. 58 

5850 

DATA 

1 4400, 

130. 27 

5860 

DATA 

21600, 

67.97 

5870 

DATA 

28800, 

42.  48 

5880 

DATA 

36000 , 

28.  32 

5890 

DATA 

99000, 

28.  32 

5900 

DATA 

9999 , 0 

Cl  2 


APPENDIX  D:  NOTATION 


A  Wetted  cross  section  of  flood  channel 

A^  Wetted  cross  section  of  dam  breach 

A  Surface  water  area  within  reservoir 
s 

A^  j  Partial  area  of  sliding  wedge 
b  Bottom  width  of  the  breach 
B  Top  width  of  the  breach 
Bp  Top  width  of  the  dam 
c^  Wave  celerity 
C  Cohesion 

Cp  Chezy's  coefficient  of  friction 
C.  Muskingum-Cunge  coefficients  (i  =  1,  5) 

* 

C  Broad-crested  weir  discharge  coefficients  (i  =  1,  2) 
0^  Integration  constant 
+ 

C~  Wave  characteristics 
d  Depth  of  the  breach 

Dg  Representative  size  of  sediment  particles 

D,.q  Median  size  of  sediment  particles 

e^  Erodibility  index 

F  Breach  Froude  number 

F  Seepage  force 
H 

g  Acceleration  due  to  gravity 

G  Weight  of  sliding  wedge 

h  Hydraulic  head 

H  Reservoir  water  level 

H  Initial  reservoir  water  level 
o 

i  Imaginary  number 
I  Inflow  into  a  channel  segment 
I  Inflow  discharge 
+ 

J  Riemann's  quasi-invariants 

K  Muskingum  parameter 

K  Erosion  proportionality  constant 
c 

K  Einstein-Brown  formula  constant 
E 

l  Length  of  breach  in  the  flow  direction 


DI 


I  Minimum  breach  horizontal  length 
s 

m  Number  of  nodes  of  each  finite  element 
M  Mass  of  eroded  material 

n  Manning's  coefficient  of  friction 
N  Shape  function  (i  =  1,  2,  3,  4) 

0  Outflow  from  a  channel  segment 

p  Soil  porosity 

q  Unit  width  discharge 

q  Bed-load  discharge  in  weight  per  unit  width 

bw 

q  Lateral  inflow 

no 

Q  Total  outflow  discharge 

Q.  Breach  outflow  discharge 

b 

Qq  Crest  overtopping  discharge 

Qp  Peak  outflow  discharge 

Qg  Sediment  discharge 

Outflow  discharge  from  spillway  a< J  powerhouse  outlet 

Hydraulic  radius 

s  Side  slope  (1V:sH) 

S  Reach  storage 

S^  Energy  gradient 

Sp  Breach  shape  factor 

S  Slope  of  the  channel 

o 

t  Time 

t.  Dam  failure  duration  time 

d 

T  Top  width  of  flood  channel 

u  Water  velocity 

V  Reservoir  water  storage  capacity 

V  Water  volume  stored  in  reservoir 

x  Distance 

Xp  Horizontal  projection  of  the  sliding  wedge 
y  Water  depth  in  flood  channel 

y^  Critical  depth 

y^  Depth  of  tailwater  section 

Z  Breach  bottom  elevation 

Initial  dam  height 

u  Weighting  factor  in  Muskingum  method 
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Coefficient  of  discharge  formula 

a?  Coefficient  of  erosion  rate  formula 

6^  Exponent  of  discharge  formula 

$2  Exponent  of  erosion  rate  formula 

Y  Specific  weight  of  water 

y  Specific  weight  of  soil 
s 

Yj  Specific  weight  of  saturated  soil 
Y0  Specific  weight  of  submerged  soil 
At  Time  step 

Ax  Length  of  channel  segment 
AZ  Depth  oc  scour 

C  Angle  between  shearing  plane  and  horizontal 
0  Angle  between  breach  sides  and  vertical 
v  Kinematic  viscosity  of  water 
t  Shear  stress 

<j>  Angle  of  repose  of  the  soil 
$  Sediment  transport  rate  function 

(g) 

X  Unknown  variable  within  a  finite  element 
Unknown  variable  at  nodal  point 
T  Inverse  of  Shields  dimensionless  shear  stress 


D3 


